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Abstract 

Gaufi cubature (multidimensional numerical integration) rules are the nat- 
ural generalisation of the ID Gaufi rules. They are optimal in the sense that 
they exactly integrate polynomials of as high a degree as possible for a par- 
ticular number of points (function evaluations). For smooth integrands, they 
are accurate, computationally efficient formulae. 

The construction of the points and weights of a Gaufi rule requires the so- 
lution of a system of moment equations. In ID, this system can be converted 
to a linear system, and a unique solution is obtained, for which the points 
lie within the region of integration, and the weights are all positive. These 
properties help ensure numerical stability, and we describe the rules as 'good'. 
In the multidimensional case, the moment equations are nonlinear algebraic 
equations, and a solution is not guaranteed to even exist, let alone be good. 
The size and degree of the system grow with the degree of the desired cuba- 
ture rule. Analytic solution generally becomes impossible as the degree of the 
polynomial equations to be solved goes beyond 4, and numerical approxima- 
tions are required. The uncertainty of the existence of solutions, coupled with 
the size and degree of the system makes the problem daunting for numerical 
methods. 

The construction of Gaufi rules for (fully symmetric) n-dimensional regions 
is easily specialised to the case of C/3, the unit sphere in 3D. Despite the 
problems described above, for degrees up to 17, good Gaufi rules for U3 have 
been constructed/discovered. 
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1 Introduction 



1.1 Multidimensional Gaufi Cubature 

Instead of directly considering the surface of the unit sphere U3, we will consider a 
more general case. For n-dimensional regions TZn, we will construct A^-point cuba- 
ture rules {xi, Wi}^^^ of the form: 

N 



/ '^(x)/(x)dx w V'w;i/(x, 



4=1 

The Xi G are called the cubature points, and the G R are their respective 
weights. We want rules that are as accurate as possible for a given number of 
points (function evaluations). Smooth integrands may be accurately approximated 
by polynomials, hence GauB rules, which exactly integrate all polynomials of as high 
a degree as possible, are a natural choice. For numerical stability, we also want rules 
with positive weights, and points which lie within TZn (these things are automatic 
in ID). Up until about 1975, the best rules known were the (Cartesian) product 
rules, which although good, are not optimal. 

The theory behind multidimensional GauB rules for 'fully symmetric' TZn is a 
natural generalisation of the well-known ID case. It originates in a foundation paper 
by Mantel and Rabinowitz (1977) ||l^. The material has been applied to the case 
of C/3 [D, although no computed rules have been pubhshed. 

We describe, in §||, the theory behind the construction of GauB cubature rules 
for fully symmetric regions TZn- The material is specialised to the case of in 
§1^, which describes implementational details and some computer programs. For 
degrees up to 17, good cubature rules for Uz have been found, and these are listed 
in Appendix ^ By comparison, the (non-optimal) product rules for C/3 are simple 
to construct, and are not terribly inefRcient. We describe them, and some programs 
for their computation, in Alternative methods for cubature, such as Monte Carlo 



and lattice methods are discussed in §4.4. In general they are poor second choices 
when GauB rules are available. 

Before progressing, we sketch an application of the use of cubature for C/3: the 
solution of the interior Dirichlet problem using a boundary integral equation. More 
significant applications are from statistics, where cubatures over many dimensions 
are required, and efficiency is critical. 



1.2 Application: 3D Interior Dirichlet Problem 

Definition 1.1 (3D Interior Dirichlet Problem) Given a smooth, hounded do- 
main G C K'^, with boundary dG, find u : G — s- R such that: 

1. u satisfies Laplace's equation within G, that is V^ii(x) — 0, Vx G G. 

2. u is known on the boundary (the Dirichlet condition). That is, there is a 
continuous Junction f : dG R such that: Vx G dG, u(x) = /(x). 

We will consider only domains G C R'^ of class G^ |l^, pp 21-22], which we 
will loosely call 'smooth'. (Their boundaries dG will be of class G^.) By G^{G), we 
mean the set of k times continuously differentiable real- valued functions defined on 
G. We will be interested in functions u contained in G(G) n G^{G). 
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The interior Dirichlet problem is a convenient model problem to work with. Its 
solution represents a potential function, that is readily related to physical observ- 
ables (e.g. electrostatic force). It is mathematically attractive, as the existence of a 
unique solution is known. Where the boundary 9G, and the boundary data / are 
simple, it may even be possible to find an analytic solution. In general, however, 
this is not possible, and it is more sensible to construct a numerical approxima- 
tion. Knowledge of the existence of a unique solution, for even quite nonsmooth 
boundaries, greatly encourages this. Techniques to (approximately) solve the inte- 
rior Dirichlet problem may be able to be used as models for the solution of more 
sophisticated boundary value PDEs. 

One method of solving the interior Dirichlet problem is to reformulate it in terms 
of a boundary integral equation on dG. This reformulation is a natural choice: it 
is involved in a constructive proof of the existence of a solution . We begin by 
defining the 'fundamental solution of Laplace's equation in 3D' as the function:]^ 



For fixed y G M'^, $(-,y) is harmonic (satisfies Laplace's equation) in \ {y}. 
Writing n(y) as the unit outward normal of dG at the point y, we construct the 
solution to the interior Dirichlet problem in terms of the fundamental solution as 
follows: 

Theorem 1.1 (Solution to the interior Dirichlet problem) For x e G, the 

double layer potential 



with continuous density (f> is a solution of the interior Dirichlet problem if (p is the 
solution of the following integral equation, for x G dG: 



A numerical approximation to the solution (f) of the BIE can be used to con- 
struct a numerical approximation to the double layer potential u, and hence the 
solution to the Dirichlet problem . Initially, we construct a cubature rule for ap- 
proximating integrals over dG. This cubature rule is used with a Galerkin technique 
to approximate the solution of the BIE for (jj. Lastly, u is approximated within G 
by application of the cubature rule to the Galerkin approximation. A cubature rule 
for the manifold dG may be constructed by pointwise projection from one for U3. 
An advantage of this is that for a particular number of desired points in the rule, a 
single rule for U3 will suffice for any manifold. A disadvantage is that although the 
rule may be appropriate for the mapping process may reduce its efficacy. An 
even point density on U3 may lose its regularity when mapped, especially if dG is 
not concave. To illustrate this, |^ uses the following region in experiments: 

1 / acos(0) sin(0) \ 

X = J cos(26l) + Jc- sm'^{20) b sin(0) sin(6l) 
^ V cos{9) J 

A typical choice is (a, 6, c) — (1,2,1.1). As c decreases towards 1, the shape be- 
comes less convex, eventually becoming like a peanut, and numerical methods lose 
accuracy. 

^ A different function is defined for a different number of dimensions, notably 2. 




</>(x) - 2 / ^(y)^l^ds(y) = -2/(x) 
JdG On(y) 
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2 Gaufi Cubature for Fully Symmetric Regions 



2.1 Motivation for Cubature 

Let TZn be an n-dimensional region contained in R"; and lu, f : TZn — >■ R- Consider 
the numerical approximation of multiple integrals of the form: 



'^(x)/(x)dx. 



Here, w is a weighting function, typically defined to contain the 'singularity' of the 
integrand. For example if TZn is K^, a typical weighting function is — e^''^' , and 
for a wide class of functions / (e.g. polynomials), the integral will exist. Commonly 
to will be unity, and TZn will be spherically symmetric about the origin. For example, 
if TZn is U3, we might have: 



/ xiX2x'^ dxi dx2 dx3 — 0. 



'Us 

Accurate approximation of multidimensional integrals is in general a computation- 
ally expensive task. The rapid growth in expense with n has been called the 'curse 



of dimensionality' (see §4.1). Fortunately, for many applications, TZn has some sym- 
metry (consider U3), and this can greatly simplify both the construction and the 
application of rules. Here, we consider the case where TZn bas some symmetry (not 
an issue in ID!), and the construction of GauB rules, which minimise computational 
expense for a desired accuracy. The philosophy is that if TZn and uj have a certain 
symmetry, then it will be natural to place equally-weighted cubature points within 
R" according to the same symmetry. This symmetry should considerably reduce 
the computations required to construct the rule. 



2.2 Full Symmetry 

One of the most natural symmetries to conceive for R" is full symmetry in Cartesian 
coordinates. A set of points is fully symmetric if any point can be reached from any 
other point by a series of orthogonal rotations about coordinate axes, and reflections 
in coordinate planes. Observe that in 3D, the set of vertices of the unit octahedron 
is fully symmetric; indeed a complete fully symmetric set of points is an orbit under 
the action of the group Gg of symmetries of the octahedron. This group has order 
48, and includes both rotations and reflections. In particular U3 is composed of 
complete sets of fully symmetric points. 

The use of symmetry for cubature rules dates back to Russian works in the 
1960s, primarily in papers by Sobolev ^ Other symmetries for 3D can 
be defined in terms of the symmetry groups of other Platonic solids. The immedi- 
ately attractive case is that of the dodecahedron/icosahedron, but this is harder to 
visualise than that of the octahedron (where the vertices all lie on the coordinate 
axes), and leads to more difficult algebra to disentangle. In any case, alternative 
symmetries have little application beyond 3D, whilst the notion of full symmetry 
generalises perfectly. The notion of full symmetry dates back to Lyness (1965) 
and the use of octahedral symmetry for cubature on U3 was first put on a clear 
foundation by Lebedev (197 6) fl^ . Shortly after this, the foundation paper of Man- 
tel and Rabinowitz (1977) ||l8[ generalised the notion to arbitrary fully symmetric 
domains, implicitly using the octahedral symmetry, although not acknowledging the 
intellectual heritage. This presentation of the theory closely follows both jl^ and a 
complementary paper by Keast and Lyness (1979) To begin, we formalise the 
notion of full symmetry. 
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Definition 2.1 (Full Symmetry between Two Points) x and y are a pair of 
fully symmetric (FS) points in R", denoted x ~ y, if y can be reached from x by 
permutations and/or sign changes of the entries ofx.. 

Observe that "~" is an equivalence relation, and thus it induces a partition of 
R" into equivalence classes. Suppose that point y has r non-zero coordinates, of 
which p are distinct, and let the jth distinct non-zero coordinate appear Ij times 
(j = 1, . . . so that h + . . . + Ip = r. Trivially, ^ p ^ r ^ n. For this y, a point 
X ^ y can then be found such that: 

■X. = {xi,Xi, ... ,Xi,X2,X2, ■■■ ,X2, ... ,Xp,Xp, ... ,Xp, 0,0, ... ,0). (1) 

/itimes /atimes iptimes 

Here we have chosen < xi ^ X2 ^ . . . ^ Xp. This x is called the generator of the 
equivalence class, and each class can be expressed uniquely in terms of a generator. 
The number of distinct points in the equivalence class containing generator x can 
be shown to be: 

2"n! 

{n-ry.hU2\...lpl' 

The number of elements in each class thus varies between 1 (for the equivalence 
class [0]) and 2"n! (in general). The exact number will be important for our work. 

Definition 2.2 (Full Symmetry for a Set of Points) A set of points TZn C R" 
is called fully symmetric (FS) if ^ TZn y ^ x imply y G TZ^, that is the set 
TZn contains only complete equivalence classes. 

Clearly R" is an FS set of points. Other important examples of FS sets of points 
(domains of integration) commonly found in the literature are the n-dimensional 
unit hypercube C„ = [—1,1]"; the n-dimensional unit sphere 

Sn = {x e R" \ xf + . . . + x'^ ^ l}; and its surface C/„. In f||, we will be interested 
in U^. 

Definition 2.3 (Full Symmetry for a Function) Given an FS domain TZn, a 
function g : TZn — > R is said to be a fully symmetric (FS) function if: for all 
x,y G TZn, X ^ y means that g{x) — g{y). 

It is not really necessary that the function be real-valued for this definition to 
make sense, but this will be sufficient. We will only consider integrals over FS regions 
involving FS weight functions, and will be approximating them using cubature rules 
on FS sets of points. Where TZn is an n-dimensional FS region, consider an integrand 
/ : TZn which is not necessarily fully symmetric, and an FS weight function 

uj : TZn which is positive over a set of positive volume: 

I[f]^ [ c.(x)/(x)rfx. 

This is to be approximated by an A^-point cubature rule {xi,i(;i}^j^, such that 
= {xii, . . . , Xin)^ G IR" and m; G R, for i = 1, . . . , A^. The rule is then: 



N 



(2) 



Recall that the degree of a polynomial in n variables is the maximum sum, of the 
exponents in any of its terms, not the maximum exponent of any one variable 
appearing in its terms. Thus 3a;f a;| + x\x3, is a polynomial in 3 (or more!) variables, 
that is of degree 7 (not 5), whilst x\ and monomials of degree 5. The 

rule is exact for a function /, if I^[f] = /[/]. If it is exact for all polynomials 
in n variables of degree up to and including k, then it is called an integration 
rule of degree (of exactness) k. Cubaturc rules of high degrees of exactness should 
accurately integrate smooth functions (which can be accurately approximated by 
polynomials), but this does not necessarily carry over to non-smooth functions. 

Definition 2.4 (Fully Symmetric Integration Rule) A cubaturc rule is 
called a fully symmetric integration rule if the evaluation points form an FS set, 
and all points in an FS equivalence class within the rule have the same weight. 

Note that this does not rc^quire / to be FS, only that the set of evaluation points 
X, are an FS set (and that the weights are constant over all members of the same 
equivalence class). An FS integration rule is then completely specified by a set of 
generators and their corresponding weights. This can greatly simplify computations 
by reducing the number of points in a rule. The degree of an FS integration rule 
can be related to properties of polynomial integrands /: 

1. If / is a monomial containing an odd power of some coordinate variable, then 



2. If / is a monomial containing only even powers of variables, then /[/] and 
[/] depend only on the exponents and not on the ordering of the variables. 

Thus, an FS integration rule which is exact for all monomials of degree up to and 
including 2m is actually a rule of degree 2m + 1, and it suffices that it be exact for 

all monomials of the form: 



where ^ < n and 1 < /c, < kj for i ^ j and fci + . . . + fc^ ^ m. (Set the monomial 

to Hi 11 = 0.) 

This rule can be written, for an appropriate set of generators X C {xj}^^, with 
elements x, each of which has an equivalence class [x] with elements y: 



xex ye[x] 
That is, to evaluate I^[f]: 

1. For each equivalence class [x], sum the function values over all elements of the 
class. 

2. Multiply the weight associated with each generator by its respective equiva- 
lence class sum, and sum these multiples. 



[/] = o. 




(3) 
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2.3 Gaufi Cubature for Fully Symmetric Regions 



Cubature rules that have a minimal number of points for a specified degree are 
called 'minimal rules'. 

Definition 2.5 (Fully Symmetric Minimal Rule) An N -point FS rule of de- 
gree 2m + 1 over an FS set TZn is a fully symmetric minimal (FSM) rule if no other 
FS rule of degree 2m + 1 over TZn exists with less than N evaluation points. 

Note that an FSM rule is not necessarily unique. In ID, FSM rules are the 
unique Gaufi rules, for which the theory is well-known. We must be careful not to 
confuse Gaufi rules with product rules (see §^ , where Gaufi-Legendre rules are used 
as basic rules in the construction of (Gartesian) product multidimensional cubature 
rules. Sometimes these rules are called Gaufi product rules. Whilst they are Gaufi 
rules, in the sense that they exactly integrate polynomials of as high a degree as 
possible (but only in one dimension), they are not minimal. We would like our rule 
to possess a couple of important properties: 

Definition 2.6 (Good Rule) An integration rule {x.i,Wi} over TZn is a good rule 
if its evaluation points lie within TZn, and its weights Wi are positive. 

The first of these conditions is familiar from ID Gaufi quadrature, whilst the 
second is new. These properties seem natural, however there is nothing in the as- 
sumptions for cubature that demands them. In ID, nai've approaches to quadrature, 
such as the (equally-spaced) Newton-Gotes family, do not preserve the positivity of 
weights for larger iV, whilst both these properties are satisfied by ID Gaufi rules. 
In higher dimensions, Gaufi rules are not necessarily good. The properties of good 
rules, in particular the latter, assist numerical stability. Whilst a rule may theoret- 
ically exactly integrate the relevant polynomials, in practice, summation of large 
terms of alternating sign tends to reduce numerical precision. In we consider the 
case where TZn is Uz, so the internal points condition will simply require points to 
lie on the surface of the unit sphere. 

We use the concept of full symmetry to define several classes of rules for the 
cubature of integrals with FS weight functions where the moment^exisi (commonly 
a; = 1). A fully symmetric cubature rule of degree 2m -\- 1 exactly integrates all 
polynomials of degree up to and including 2m. A fully symmetric minimal (FSM) 
rule does so using a minimal number of cubature points. A fully symmetric good 
(FSG) rule does so where all weights are positive and all points lie within TZn- A fully 
symmetric minimal good (FSMG) rule is both minimal and good. A fully symmetric 
good minimal (FSGM) rule is a good rule that is minimal, that is, although FS rules 
on fewer points may exist, none are good. 

We seek firstly FSMG rules, and if there are none of these, FSGM rules, which 
always exist. At worst, FSGM rules are product rules (see §^), which (generally) 
require [2{2m-\- 1) — 1]" — {Am-\- 1)" points to be of degree 2m + 1. For low- 
dimensional applications (small n), this may not be terribly inefficient, e.g. product 
rules for actually require 2(m + 1) points. Apart from such considerations, a 
rule that is minimal (or almost so) might be almost good in the sense that negative 
weights are very small, or that points are only just outside TZm or that failing these 
conditions, errors in integrating polynomials are minor. The quest for good rules 
may be an arduous search through these almost good rules. 

^ That is, the integrals of appropriate polynomials over 7?.n- 
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2.4 Conditions for Gaui3 Cubature in 3D 



The procedure for constructing Gaufi rules becomes more complicated as the di- 
mension n increases. Here, we fully develop the 3D case, specialising this to the 
case of U3 in §||. We will set up a system of (moment) equations to express the fact 
that an FS rule for TZ^ will exactly integrate all polynomials of up to a specified 
degree, without presupposing the number or distribution (beyond being an FS set) 
of points, or the sign of weights. 

To begin, we group generators into types, depending on the number of zeros and 
repeated elements in the entries of their equivalence classes. There will always be 
the (not very interesting) class [0] , of unit multiplicity, and in general e other types 
of classes, for a total of e + 1 types of classes. For small n, this e is generally small, 
and may be found as the solution to the following problem:^ 

Given a positive integer n, e + 1 is the total number of strings of length up to 
n, with p distinct entries taken from the positive integers 1, . . . ,p, of the form: 

2,2,^. . ,2 , . . . , p,p, . ^. . ,p P^n, 
iitimes /2times /ptimes 

such that h ^ . . . ^ Ip and li + . . . + Ip ^ n. 

Answers to this problem can be found by counting the strings, and this is 
implemented in C as findec.c (Appendix Output from this program (for 
n = 1, . . . , 100) is presented in Appendix Not surprisingly, e grows rapidly with 
n. (This is just for curiosity purposes; will only apply the case n = 3.) 

For the case n = 3, there are 7 types of classes of points, listed in Table |l|. There 
could be a generator at the origin, so that's one type, called type [0]. Generators on 
a coordinate axis are of a second type, called type [1], in which there are 6 members 
of each equivalence class. Generators of the form {(3,(3,0), in a class of size 12, are 
of type [1,1], etc.^ More generally, the generator in (j^) is of type [li,l2, ■ ■ ■ ,lp], and 
the complete set of types of generators in 3D is included in Table This notation 
(from 1^ ) , is simplified for the case n = 3 in , where there are Ki generators of 
each type, for i = 0, . . . , 6, and we shall use the latter notation. We shall refer to 
our rule as having structure {Ki}^^^ (an ordered set), usually just written {Ki}. 



Class 


Class 


Number of Generators 








Class 


Number 


Type 


([§ and 11) 


Names of Generators and Weights 


Size 





[0] 


K[0] = Ko 


(0,0,0), 





if Ko = l 


1 


1 


[1] 


K[l] = Ki 


(a.,0,0). 




i = l,...,Ki 


6 


2 


[2] 


K[2] = K2 




h 


l=l,...,K2 


12 


3 


[1,1] 


K[l, 1] = K^ 




Ci 


i = l,...,K3 


24 


4 


[3] 


X[3] = Ki 




di 


i = l,...,Ki 


8 


5 


[2,1] 


K[2, 1] = K^ 


(Ci , Ci , ) , 


ei 


l = l,...,K5 


24 


6 


[1,1,1] 


K[l,l,l]=Ke 


{9i,fJ.t, Xi), 




i = l,...,Ke 


48 



Table 1: Nomenclature of generators and weights for Gaufi cubature in 3D. 



^ The framing of this problem is the basis for the higher-dimensional analysis in 

* Each of these generator types can be thought of in terms of geometrical arrangements of 

points on t he s urface of a unit sphere, projected from the vertices, edges and faces of the unit 

octahedron |14[. 
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Given a structure {Ki}, there will be a total of A'^ points in our rule, given by: 



N^Ko + 6Ki + 12K2 + 24:Ks + 8K4, + 24X5 + ASKq. (4) 

To construct GauB rules, we proceed without presupposing the structure, instead 
attempting to find a structure such that the conditions for GauB cubature are 
satisfied, and N is minimised. Taking m (such that the desired degree is 2m + 1), 
and a rule structure {Ki}, we write down a system of (moment) equations involving 
an appropriate set of generators, based on the requirement that a GauB rule exactly 
integrates all n-variable polynomials of each degree up to 2m. It is in fact a sufficient 
requirement that we integrate exactly (n-variable) monomials of these degrees. The 
system of moment equations is a system of nonlinear algebraic equations. (This is 
also true in ID, although clever artifice allows us to reduce its solution to that of a 
linear system.) 

For FS TZs, the variables are listed in Table |l|, and the system (^) expands 
to that presented in Figure 0. We will call this system (*), in accordance with 
pp 410-411], where it first appears explicitly. System (*) splits naturally into 
three subsystems, defined by the number of non-zero indices ki in the monomial 
ccj'^^ ccg*^^ Xg*^^ that we wish to integrate exactly. Recall that the {Ki} are non-negative 
integers, and the other variables are real. Also, although Kq G {0, 1}, we include it 
in a sum for consistency, and apply the convention that o = if Kq = 0. If 

m < 3, subsystem III is ignored, and if m < 2, subsystem II is also ignored. 

Row-wise examination of the variables in Table shows that, given {Ki}, there 
are a total of: 

v^Ko + 2Ki + 2K2 + 3K3 + 2Ki + 3if5 + 4^6 

variables in (*). To determine the number of equations, let r be a positive integer 
and, for v — \, ... ,n, lei pu{r) be the number of solutions in positive integers ki of: 

ko + ki + . . . + ky = r 1 ^ fco ^ fci ^ . . . ^ fci,. 

For z/ = 0, we say there's one solution if r = 1, and none otherwise. Pu{r) is the 
dimension of the space spanned by: 



x\^^ ■ ■ ■ xl^" \ ki ^ k2 ^ ks ^ ■ ■ ■ ^ ky ^ and ki + k2 + . . . + k^ 



1^ = 



This can be efficiently computed by: 

1 r = 
else 

Pu{r) = S = 1 

r + l(niod2) u ^ 2 

p^-i{r - 1) + py{r ~ ly) else. 

For each r = 1, . . . ,m, there are po{r) + pi{r) equations in subsystem I, P2ir) in II 
and P3{r) in III. Hence, for a rule of degree 2m -\- 1, the total number of equations 
involved is X^I^i Siy=o ^ ~ 1; ■ ■ • ; 20, these numbers are listed in Table 
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Subsystem I: 



Subsystem II: 



Subsystem III: 



Ka Kt K2 K3 Ki K5 Ke 

^o + 6^a, + 12^6, + 24^c, + 8^d, + 24^e, + 48^/, 

i=l 2—1 1—1 i—1 i—1 2—1 2—1 

Ki K2 K3 K4 



1=1 
K5 



i=l 



4=1 



j = 1, . . . ,TO 



K2 



K3 



= 4E ^./^^^-^-^'^ + 4E c.(7f^5f + 7^-5?^) + 8E '^.^ 



2j"+2fe 



i=l 



i=l 



i=l 

1 s$ j A: j + k^2, 



, m 



= 8 E + 8 E e. (C^^+^\f ' + + Cf '^^^^^ ) 

= 1 i=l 



-8E /^(^- Af + 0f Af + 0f Af + e^^'^' + €^^>^ + Mf Af ) 

\ i j + k + I ^ 3, . . . ,m 



i=l 



Figure 1: The system (*) of moment equations used to determine GauB cubature 
rules for regions 7^3 Q (after pp 410-411]). 



10 



r 


Po('') 




P2{r) 




J27=iJ2t=oPAr) 


1 


1 


1 








2 


2 





1 


1 





4 


3 





1 


1 


1 


7 


4 





1 


2 


1 


11 


5 





1 


2 


2 


16 


6 





1 


3 


3 


23 


7 





1 


3 


4 


31 


8 





1 


4 


5 


41 


9 





1 


4 


7 


53 


10 





1 


5 


8 


67 


11 





1 


5 


10 


83 


12 





1 


6 


12 


102 


13 





1 


6 


14 


123 


14 





1 


7 


16 


147 


15 





1 


7 


19 


174 


16 





1 


8 


21 


204 


17 





1 


8 


24 


237 


18 





1 


9 


27 


274 


19 





1 


9 


30 


314 


20 





1 


10 


33 


358 



Table 2: The total number of equations involved in system (*), for m = 1, . . . , 20. 



In general, (*) may not have a unique solution, and indeed, may have no solu- 
tions. The best analytic technique currently available for the solution of systems of 
multivariable polynomial equations is the use of Groebner bases, see for example 
1^. (The computer algebra package Mathematica uses this method.) Symbolic 
computations are expensive, however, and will fail to yield answers for polynomial 
equations of degree greater than 4. We are quickly led to numerical techniques! For- 
tunately, experience with (*) shows that solutions commonly do exist; we progress 
with this as a hope. 

For a given degree 2m + 1, there are many different structures {Ki} that will 
lead to a system (*). To deduce a possible structure, we would like to choose {Ki} 
to minimise iV, such that (*) is consistent. This is an optimisation problem, with 
constraints that will ensure the loosest possible consistency of (*). It turns out pl| , 

that the optimisation problem has linear constraints. Clearly the cost function 
(m is linear, and the solutions {Ki} are integers. The optimisation problem is thus 
a linear integer programming problem, with linear cost function. This is routine to 
solve, and many examples of solutions are provided by Rabinowitz et al. in (for 
2D) and 0, |l| (for 3D). 
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We will call the constraints on (*) consistency conditions; their number will 
grow with dimension n. Establishment of the consistency conditions is in general a 
formidable problem, and can only be done manually for two or three dimensions.^ 
For 3D, it initially appears that there could be as many as 2^ — 1 = 127 constraints, 
and in general, where Tin has e + 1 types of FS equivalence classes, there might 
be 2^"*"^ — 1 constraints. Fortunately, symmetries in (*) reduce this number to a 
manageable level. For n ~ 3, this reduction can be done manually,^ and leads to 
the system of 13 constraints presented in Figure ^ The data presented in Table ^ 
can be used to convert these constraints to a linear vector inequality. 

Once the integer programming problem has been solved for a candidate structure 
{Ki}, we write the specific version of (*) for this structure, and attempt to solve it 
for generators and weights. The moments on the LHS of (*) are available analytically 
for a wide range of FS TZs (in particular U3). There are a number of subtleties in 
this process: 

1. Whilst the IPP will always have a solution, this is not guaranteed to be unique. 
We will in general have a (small) number of possible rule structures to choose 
from, and we order them lexically. 

2. Unfortunately, (*) is not guaranteed to be consistent for any particular rule 
structure deduced from the conditions (remember that they were chosen as 
the loosest possible), so we may have to try a number of possible structures 
before we succeed. 

3. Assuming that there is a solution for a particular rule structure, this is not 
guaranteed to be unique; typically there will be either a (small) finite, or an 
infinite number of solutions to {*). 

4. Rules deduced from solution of (*) may not be good. They may still be ac- 
ceptable, in cases where the points are not in but they are 'close' to it; 
or where some weights are negative, but small in magnitude. For purity, we 
reject such solutions. 

5. If there are no rules for minimal N, or those that do exist are not good, then 
we search for (FSGM) rules of a larger N. To do this, we restart the IPP 
with an added constraint that N be larger than the rejected solution. The 
new IPP will be solvable, and we continue with the solution of (*), iterating 
this procedure until we arrive at a good rule. This process is guaranteed to 
terminate, the worst possible case being that we actually construct a product 
rule (which always exist). 

6. The numerical solution of (*) may be computationally intractable. It may 
be difficult to tell whether an algorithm is not converging because a solution 
does not exist, or because the system is so severely nonlinear that the software 
cannot handle it. 

Despite these problems, solutions to (*) have been computed for several TZn, 
and various weighting functions. Primary results for TZ2 and TZs are contained in 
[ p7| , |l8|, . Some results for U3 (see §^) , are presented in Appendix ^. 



^ Although a framework for the automatic construction of consistency conditions in any number 
of dimensions is presented in [p^, it seems to he essential to use a machine to do this. 
^ See [[|, pp 394-398], and extra details in [0. 
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Figure 2: The constraints for system (*), expressed in terms of the Pu{r) (after 
p 398], in which it is caUed system C). We apply the convention that J2i=p'^i — 
if q < p. 
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3 Gaufi Cubature for 



We specialise the theory presented in §g to the case where TZn is f/3, and the weight- 
ing function is unity. As the 'volume' of C/3 is finite, integrals over it of any bounded 
function will be defined. We describe the spherical harmonic polynomials, which our 
cubature rule should integrate exactly, and comment that this is equivalent to our 
previous requirement that our rules integrate polynomials exactly. The construction 
of GauB cubature rules follows, with simplifying assumptions that are obtained by 
considering the region of integration. The material presented is related to that of 
Keast 0, |. 

3.1 Spherical Harmonics 

Spherical harmonic polynomials are the natural generalisation of orthogonal poly- 
nomials to the surface of the sphere. They satisfy Laplace's equation j?], appearing 
in particular when separation of variables is used to solve the interior Dirichlet prob- 
lem. They are orthogonal, and hence are a useful basis for approximating functions 
on U^. From p9[ |, we take the following definitions and properties. Any point x G C/3 
can be uniquely characterised by the ordered pair {9, </>) of coordinates of longitude 
and colatitude (j). The spherical harmonic polynomials Yim{0,(f>), where m ^ 
and l,m £ E, are defined in terms of the associated Legendre polynomials Pf^ as 
follows: 



They are orthogonal functions, normalised such that the integral over f/3 of a prod- 
uct YimYiim' is unity only if ^ = and m = m' . Using * to denote complex conju- 
gation: 



we can relate any spherical harmonic to an associated Legendre polynomial P™ 
with TO ^ 0. Table |^ lists some of the simplest spherical harmonics. As I and m 
increase, the degree of the trigonometric functions increases. Note that for to = 0, 
they are purely real. 





Using the relation 
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Table 3: The first few spherical harmonics and the associated Legendre polynomials 
(after |l|, p 195]). 



Any function P : C/3 M can be expanded in terms of spherical harmonics: 

00 00 



1=0 m=0 



Gaufi quadrature rules in ID exactly integrate monomials of as high a degree as 
possible using a fixed number of quadrature points. In the case of the cubature of 
integrals defined on C/3, the natural generalisation is that the rules must exactly 
integrate as many spherical harmonics as possible.^ As in the ID case, the cubature 
rule may not be optimal for any particular integrand, but should work well for 
smooth functions that are well approximated by (sums of) spherical harmonics. It 
turns out that exact integration of the spherical harmonics over is equivalent to 
applying the work in §^ with restrictions on the placement of points. We do not 



have to consider U3 as a special case, and continue on from §2.4 



3.2 Application of the Gaufi Cubature 

For C/3, the points of a good cubature rule (Table |l|) must lie on C/3, and so have some 
further constraints added to them. In particular Kq = and Ki, K2; K4 G {0, 1}. 
Also, if = 1 then ai = 1, if = 1 then /3i = l/\/2, and if = 1 then 
ei = l/\/3- Furthermore: 



In general for FS regions TZ^, there are a total of 13 consistency conditions (Figure 
2), but for C/3 these simplify drastically. There are only 4, and these listed in Figure 
3. For m = 1, . . . , 20, the appropriate right hand sides of the constraint equations 
in Figure ^ are presented in Table ^. 

^ The ID analogue of U3 is the unit circle U2, for which the orthogonal polynomials are e*"*^ 
and the points of Gaufi rules are equally spaced. 
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Figure 3: Constraints for the integer programming problem, where TZn is U3 
(c.f. Figure H). In addition, Kq = 0, and Ki, K2, K4 ^ 1. Again, the convention 
ai = if (7 < p, is followed. Values for the right hand sides are tabulated for 
some choices of r in Table 0. 
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Table 4: The four elements of the columns of the right hand side in Figure ||, for 
various m, as generated by appropriate sums of the Pj/(r). 
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Given m, for which we wish to construct a rule of degree 2m + 1, we firstly set 
up the integer programming problem to be solved for a rule structure {Ki}. As in 
Table ^, the sum of the Pv{r) over v = 0, . . . , 3, and r — 1, . . . , m gives the total 
number of equations in (*) for U^. 

Lebedev published an important paper on Gaufi cubature for just prior 
to that of the more general one by Mantel and Rabinowitz ||l^. Lebedev makes 
some astute choices in the presupposition of rule structures, which sometimes result 
in optimal choices. The algebra is simplified by enforcing Ki = K2 = = 1, and 
sometimes also K2 ~ 0; and choosing K^^Ki^, and Kc, such that the number of 
unknowns in (*) is equal to the number of equations. The rules generated were a 
great improvement on the previously completely unsystematised collection of known 
rules (e.g. see Stroud @), but we are interested in the more general case. 
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Table 5: Basic solutions to the IPP, with the integer constraint removed. Any 
integer solutions will have N at least equal to that displayed here. 
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3.2.1 Solution of the Integer Programming Problem for Uj, 

For our chosen m, we set up and solve the integer programming problem (including 
the constraints Ki, K2, K4 ^ 1, and that Kq is set to 0), for the minimisation of TV, 
the number of points in the rule: 

N = 6Ki + 12K2 + 24X3 + 8X4 + 24if5 + iSKe. 

We require all the solutions to the IPP which have N at least equal to the (integer) 
minimum, and less than some (user-specified) small multiple of this minimum. We 
firstly solve the programming problem without integer constraints. This yields a 
lower bound for N, used when searching for all the integer solutions. This initial 
problem is solved using the Mathematica program IPPBasicSoln.M (Appendix 

Results (lower bounds for N, and associated 'pseudostructures' {Ki}), for m = 
1, . . . , 20 are presented in Table ||. 

To find all the solutions of the IPP, we use a simple exhaustive search. Whilst 
this is crude, in this case it is reasonably efficient. For U^^ there are only six variables 
to search through, and Ki,K2, K4 G {0, 1}. We use the lower bound on N, and as an 
upper bound we use a small multiple (say 1.5) of the lower bound. This second phase 
is implemented as a C program ipp . c (Appendix]^). This program takes the data in 
Table ||, and exhaustively searches for all integer solutions, subject to (reasonable) 
bounds K^jK^jKq < 20, and an upper bound on A^" set empirically so that we only 
collect about the first 100 solutions. The output from ipp. c is (manually) sorted, so 
as to order the solutions firstly with increasing N , and then lexically. The solutions 
corresponding to the first five integer minima, for m = 1, . . . , 10 are presented in 
Table |. 

These minima agree with those published by Keast ||, p 155 and pp 166-167]. 
(This paper presents all structures for the first 5 consecutive minima in A^, for 
TO = 1, . . . , 9, corresponding to degrees 3,5,..., 17.) Keast claims to have obtained 
FSMG rules from the structures, but does not actually describe them, but we do. 
(Keast 's paper also correctly identifies an error in the work of Lebedev ||lj, p 15], 
but fails to note that the Lebedev's paper is less general, so that results cannot be 
compared directly.) 

Following the notation of ||l^, p 400], the rules are assigned names. Say we have 
a rule of degree 2to -|- 1 for TZn with weighting function h(r). Let this rule be found 
from the jth instance (ordered lexically) of the ith consecutive minima of the IPP, 
and have structure {Kq, Ki, . . . , K^}, and total number of points A^. We will name 
this rule: 

7^f/'') : (2m + l)-i.j{Ko, ifi, . . . , K,)-N. 

Our rules for U3, with h = 1 can be labelled as: 

C/3 : i2m + iyi.j{Ki,...,KeyN. 

When (*) has multiple solutions, the labelling is not unique, but this is not a 
problem. This nomenclature provides a basis for comparison with (published) rules 
derived by other means. 
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Table 6: All the solutions to the IPP, for m — 1, . . . , 10, corresponding to the first 
five minima. (Some of these are already available in Table ^.) v is the number of 
variables in (*), for the particular structure. 
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3.2.2 Moments for U3 

Solution of {*) for U3 requires knowledge of the moments /[x^-'iy^^^z^^^], for all 
combinations of integers ^ ji ^ j2 ^ js ^ such that ^ ji + + js ^ m. 
(The moment is zero if any of the exponents are odd, and it is invariant under 
permutation of exponents.) Explicitly, we require analytical evaluation of integrals 
of the form: 

JUa 

Clearly I[x'^y^z^] — An. Using coordinates 6*, the longitude, and the latitude (not 
the co-latitude), the integral separates p 33]: 

I[x^'^y^'^z^^^]^ I cos^^'{e)sm^^^0) de / cos2Ji+2j2+i(0) sin^^'^ (</)) d0. 

These moments are calculated by USMoments .M, a Mathematica program (Ap- 
pendix^). Table [t] lists them, for m ^ 10, and this allows us to establish (*) for 
Gaufi rules of degrees 3, 5, . . . , 21. 

3.2.3 Solution of the System of Moment Equations (*) 

Having constructed tables of moments and possible structures for various to, we 
may consider solution of (the pared-down version of) (*) for U3. Attempting to 
use Mathematica to do this succeeded for to up to 5 (sometimes with a little 
human assistance in making substitutions). Beyond to = 5, the high degree of the 
polynomials involved means that (*) in general has only transcendental solutions. 
Approximate solution of (*) is attempted using numerical software in MATLAB (and 
C). 

We use the matlab routine f solve, which numerically approximates the solu- 
tion of a system of equations. This routine requires a user-specified function that 
evaluates the system (*). Initially, this was done as a MATLAB m-file. As the size of 
(*) increases rapidly with to, the large number of expensive function evaluations 
involved made this a slow procedure. Instead, momenteq. c, a MATLAB mex-file was 
written (Appendix ^), which yielded a speedup of two orders of magnitude. As 
numerous experiments are required to find a satisfactory structure, this speedup 
is important. To aid debugging, a program writestar.c (Appendix takes in- 
puts of TO and {Ki}, and outputs a M^TjX file containing (*) as a set of displayed 
equations. 

A MATLAB driver program, cubature.m (Appendix ^) is used to try various 
structures and choices of to. For each to, cubature.m is run until either f solve 
solves (*), or the user gives up in disgust. Successful results for to = 1, ... ,8 are 
presented in Appendix 
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j2 


is 




ii 


h 


is 










1 


47r/3 





4 


6 


47r/12597 








2 


47r/5 





5 


5 


127r/46189 








3 


47r/7 


1 


1 


1 


47r/105 








4 


47r/9 


1 


1 


2 


47r/315 








5 


Air/ll 


1 


1 


3 


47r/693 








6 


47r/13 


1 


1 


4 


47r/1287 








7 


47r/15 


1 


1 


5 


47r/2145 








8 


47r/17 


1 


1 


6 


47r/3315 








9 


47r/19 


1 


1 


7 


477/4845 








10 


47r/21 


1 


1 


8 


47r/6783 





1 


1 


47r/15 


1 


2 


2 


47r/1155 





1 


2 


47r/35 


1 


2 


3 


47r/3003 





1 


3 


47r/63 


1 


2 


4 


47r/6435 





1 


4 


47r/99 


1 


2 


5 


47r/12155 





1 


5 


47r/143 


1 


2 


6 


47r/20995 





1 


6 


47r/195 


1 


2 


7 


47r/33915 





1 


7 


47r/255 


1 


3 


3 


47r/9009 





1 


8 


47r/323 


1 


3 


4 


47r/21879 





1 


9 


47r/399 


1 


3 


5 


47r/46189 





2 


2 


47r/105 


1 


3 


6 


47r/88179 





2 


3 


47r/231 


1 


4 


4 


287r/415701 





2 


4 


47r/429 


1 


4 


5 


47r/138567 





2 


5 


47r/715 


2 


2 


2 


477/5005 





2 


6 


47r/1105 


2 


2 


3 


47r/15015 





2 


7 


47r/1615 


2 


2 


4 


47r/36465 





2 


8 


47r/2261 


2 


2 


5 


1277/230945 





3 


3 


207r/3003 


2 


2 


6 


47r/ 146965 





3 


4 


47r/1287 


2 


3 


3 


47r/51051 





3 


5 


47r/2431 


2 


3 


4 


477/138567 





3 


6 


47r/4199 


2 


3 


5 


47r/323323 





3 


7 


47r/6783 


2 


4 


4 


47r/415701 





4 


4 


287r/21879 


3 


3 


3 


207r/969969 





4 


5 


287r/46189 


3 


3 


4 


207r/2909907 



Table 7: Moments of the first few monomials over C/3. 
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4 Product Rules for 



Product rules are an important approach to multidimensional cubature, and they 
are the main direct competitor with Gaufi rules. Here we describe their construc- 
tion for U3, and compare their efficiency with the Gaufi rules discussed in §H The 
material is largely abstracted from chapter 2 of Stroud pGj . 

4.1 Introduction 

To introduce the concept of product rules, consider the case where TZn is C3 (the 
unit cube), and we have a unit weighting function. Let {xi,Wi}^^ be the m-point 
ID Gaufi-Legendre rule of degree 2m ~ 1, which exactly integrates all of: 
f-i 

x-^dx j = 0, . . . , 2m — 1. 

-1 

For C3, we wish to integrate exactly: 

rj^Xj^Xg^ dxidx2dx3. 

-1 j-i j-i 

Writing this as an iterated integral allows us to construct a product rule, which 
exactly integrates all of: 

1 /.I pi 



xf dxi / x^^ dx2 I xl" dx3 ji, j2, ja = 0, . . . ,2m - 1. 

-1 J -1 J-i 

Let i ~ (ji,«2,i3) for all 1 ^ zi,Z2,«3 ^ m. We construct an m'^-point rule {xi,^i} 
for C3 via: 

X; = (xij , , 2^43 ) and Ai = Wi-^Wi^wi^. 
The set of points is a Cartesian product: 

{Xi} = {Xi^} X {Xi^} X {Xi^}. 

This product rule will integrate exactly all monomials: x 2:2^ , for ji,j2,j3 = 
0, . . . , 2to — 1. The highest degree monomial that it will integrate exactly is 
xl'^^^xf^^^xf"^^ , of degree 3(2m-l). It is not however, a rule of degree 3(2m-l) 
as it does not exactly integrate all polynomials of this degree, for instance, it does 
not exactly integrate x^"^™" . It is a Gaufi rule in that it exactly integrates all 
polynomials of up to a certain degree, but only in one variable. However, as the ID 
Gaufi rules are good, this product rule is also good. This is a general property of 
product rules. 

The rule requires m"^ points. In general, a product rule on TZm of degree 2m — 1 
in each of the n variables, will require N = m" points. This exponential growth 
in N was called the 'curse of dimensionality' by authors in the 1960s, when it was 
believed that there was no escape from it. 

For small n product formulas are very useful. For example if one 
wanted a subroutine that used a fixed 1000-point formula for a wide 
class of integrands for the 3-cube, w{x,y,z) = 1, we believe that one 
could do no better than the product of three copies of the 10-point 
Gaufi-Legendre formula. This formula has degree 19 and there are, in 
fact, no nineteenth-degree formulas known for the 3-cube using fewer 
than 1000 points. (. . . a lower bound for the number of points in such a 
formula is 221.) 

Stroud (1971), g p 25]. 
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n 


m 


N = to" 


TO 


iV = to" 


TO 


N = to" 


2 


31 


961 


999 


998001 


31622 


999950884 


3 


9 


729 


99 


970299 


999 


997002999 


4 


5 


625 


31 


923521 


177 


981506241 


5 


3 


243 


15 


759375 


63 


992436543 


6 


3 


729 


9 


531441 


31 


887503681 


7 


2 


128 


7 


823543 


19 


893871739 


8 


2 


256 


5 


390625 


13 


815730721 


9 


2 


512 


4 


262144 


9 


387420489 


10 


1 


1 


3 


59049 


7 


282475249 


11 


1 


1 


3 


177147 


6 


362797056 


12 


1 


1 


3 


531441 


5 


244140625 


13 


1 


1 


2 


8192 


4 


67108864 


14 


1 


1 


2 


16384 


4 


268435456 


15 


1 


1 


2 


32768 


3 


14348907 


16 


1 


1 


2 


65536 


3 


43046721 


17 


1 


1 


2 


131072 


3 


129140163 


18 


1 


1 


2 


262144 


3 


387420489 


19 


1 


1 


2 


524288 


2 


524288 


20 


1 


1 


1 


1 


2 


1048576 



Table 8: Limits on to for various n using product rules requiring N = to" evaluation 
points, if ceilings of < 10'^, 10® and 10^ are enforced. (Modelled after table 2.1 in 
@P24].) 

In fact shows that in theory, an FSM rule for C3 can be constructed using 
345 points, although there is no explicit calculation of such. It seems likely that an 
FSGM rule on less than 400 points exists. This book illustrates the limitations 
on the usefulness of product rules in terms of the relationship = to". Table ^ 
illustrates some upper limits on to for various n if a ceiling of function evaluations 
is enforced. 



25 



4.2 Product Rules for U3 

Any point in can be uniquely characterised by the longitude 9, and the co- 
latitude (f). This simplifies the construction of product rules, as we can express such 
a rule as the product of rules that integrate over 9 and respectively. The following 
construction is abstracted from ||2^, pp 34-35 and 40-41]. 

For k — 1,2, let {yk,i, ^fc,i} be the points and weights in the m-point ID Gaufi- 
Jacobi rules: 

/I rn 
1 i=l 

For k — 1, this is the Gaul3-Chebyshev rule of the first kind {yi^i, Ai^i}: 

/I m 

For i = 1, . . . , m, these are [|l3[ p 114], given by the following formula (note that 
the weights Ai^i are constant): 

,(2i-l)7r, ^ ^ , 

yi.i = cos( ) and An = t: m. 

2m 

For fc = 2, the formula is a GauB-Legendre rule {2/2,1, ^2.1}- 



/I m 
f{y2)dy2 « VA2,,/(y2,,). 
-1 



This is not expressible in a simple closed form, but calculation is routine and efficient 
(e.g. the implementation in gauss. m in Appendix ^% . 

Now let i = {11,12), for \ ^ 11,12 ^ m, and define v-^ i, fi 2, 1^1,3 by: 

Via = ±(i-y2,»2) 

A 2rn^-point product rule {xi, B{\ of degree 2ra — 1 for [/a is then given by: 

Xi = (i^i,i, 1^1,2, 1'i.a)^ and Si = Ai^i^A2a2- 

Some substitutions and relabelling in the construction shows that it is actually 
simple, and we write the computations as an algorithm. 
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Algorithm to Construct Product Rules for U3 

1. Given to, create a matrix i = {ii, 12), where 1 ^ ii, 12 ^ to. This matrix labels 
indices of the points in the rule. 

2. Compute the TO-point Gaufi-Legendre rule {yi^,Ai^}. 

3. Set 



cos( 

±yi2 



(2ii-l)7r 
2m 



and 



m 



as the 2to^ evaluation points, and the corresponding weights of a product rule 
of degree 2to — 1. 

Implementation as two matlab files is presented in Appendix ^ Specific illus- 
tration of the points and weights of the rules generated is rather pointless, but a 
table of errors for integrating the appropriate polynomials generated by the function 
uSprod.m shows it to work perfectly (the errors are of order machine precision). 
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4.3 Comparison with Gaui3 Rules for Us 



For the the case of U3, the curse of dimensionahty for product rules is not the 
terrible scourge that it might have been. Table ^ compares the number of points 
required with varying degree for FSGM (or even FSMG) rules, as compared with 
the product rules. In many applications, we may have no reason to want rules of 
degree more than about 10. Observe that there is only a small factor of inefhciency 
in using the product rules, not the many orders of magnitude that appear when the 
dimension is higher. 





N 


Ratio 


Degree 


FSMG 


Product 


% 


3 


6 


8 


75 


5 


14 


18 


78 


7 


26 


32 


81 


9 


38 


50 


76 


11 


50 


72 


69 


13 


78 


98 


80 


15 


86 


128 


67 


17 


110 


162 


68 



Table 9: Comparison of the number of points required for the FSGM (or FSMG) 
Gaufi rules and the product rules for U3. The last column is the ratio of the number 
of points (function evaluations) required by the Gaufi rules relative to the product 
rules. 



4.4 Alternative Philosophies 

Gaufi rules (and their ancestors the equally-spaced formulae) are based on exploiting 
the analytical properties of smooth integrands. They are optimal in the sense that 
in general they will be the best choice for approximating integrals involving smooth 



integrands. As mentioned in §2.2, the optimal multidimensional cubature can only 
be sensibly considered for regions with some symmetry, and we dealt with the case 
of full symmetry. 

To deal with non-smooth integrands, possibly even random distributions, and 
with non-symmetric, possibly even disconnected domains of integration, alternative 
philosophies are usually more relevant. Textbooks on numerical integration com- 
monly contain many pages describing minute implementational details to further 
refine the theory for optimal methods; and then devote a similar amount of space 
to describing real alternatives (e.g. |6[ ^). The main alternative technique is the 
Monte Carlo method, but a more recent idea is the 'lattice' method. 



4.4.1 Monte Carlo 'Simulation' Methods 

These techniques are based on averaging function values at a random selection 
of points within the region, and they are particularly appropriate for oddly-shaped 
regions and non-smooth integrands. They are widely used in statistical applications, 
where integrals of high dimension must be approximated. Performance is often about 
N^^ (the error incurred using N points should be proportional to N^^). Volumes 
have been written about them (e.g. see |2^, chapter 6]). 
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4.4.2 Lattice Methods 

These methods generahse the idea of placement of equally-spaced points in ID (with 
weights selected according to some generally simple formula, expressible in closed 
form). The idea is to catch as representative a sample of function values as possible. 
In ID, this leads to the Newton-Cotes family of rules. For smooth integrands, these 
rules increase in accuracy algebraically with their number of evaluation points, 
although they are not always good. (Whilst this may be adequate, it is still inferior 
to the exponential accuracy of the Gaufi rules.) 

The problem with attempting to generalise this to the case of several dimensions 
is that the notion of 'equal spacing' of evaluation points becomes less well-defined. 
Placement of equally-spaced points on U3 is equivalent to maximally covering it 
with nonintersecting equal circles, a problem thought to be intractable. The 'best' 
that can be done involves heuristic algorithms, and lots of computer time Q. 

Nevertheless, it may be reasonable to try to approximate the equally-spaced 
placement of points within our region. Recent research involving Sloan and Lyness 
[ p6| , |2^, has achieved this using geometric construction techniques, called 'lattice 
methods'. For periodic functions on [0,1]^, lattice methods generalise the trape- 
zoidal rule, preserving the order of the error as N~'^. For U3, a placement called 
'spherical i-designs' is used 



4.4.3 Comparison vi^ith Gaufi Cubature 

For FS regions and weight functions, where the expected integrands are smooth, 
these alternatives are a poor second choice in comparison with Gaufi rules. The 
Monte Carlo methods will only converge as iV~^, and the theory for the lattice 
methods is not very general, results only being available for one type of region at a 
time. 

It cannot be overemphasised that where Gaufi cubature is available, it should be 
used, particularly as the dimensionality increases. For evaluation of integrals over 
2D manifolds, Gaufi cubature is applicable. 
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A Number of Equivalence Classes with Dimension 



Results from running f indec . c (Appendix ^) for n = 1, . . . , 100 are listed in Table 
p^ , and graphically presented in Figure ^. For n = 100, the (optimised) program 
requires about 1000 CPU minutes on a SPARC- 10 workstation. Data has been 
manually checked for n = 1, . . . , 10. 



n 


e + 1 


n 


e + 1 


n 


e + 1 


n 


e + 1 


n 


e + 1 


1 


2 


21 


3506 


41 


259891 


61 


7760854 


81 


141227966 


2 


4 


22 


4508 


42 


313065 


62 


9061010 


82 


161734221 


3 


7 


23 


5763 


43 


376326 


63 


10566509 


83 


185072690 


4 


12 


24 


7338 


44 


451501 


64 


12308139 


84 


211616350 


5 


19 


25 


9296 


45 


540635 


65 


14320697 


85 


241783707 


6 


30 


26 


11732 


46 


646193 


66 


16644217 


86 


276046669 


7 


45 


27 


14742 


47 


770947 


67 


19323906 


87 


314934342 


Q 
O 


D ( 


zo 


io40U 




y lozzU 


Do 


00/1 1 1 ^1/1 1 

zz41 1041 


OO 


ooyU4z4ol 


9 


97 


29 


23025 


49 


1091745 


69 


25965986 


89 


409038376 


10 


139 


30 


28629 


50 


1295971 


70 


30053954 


90 


465672549 


11 


195 


31 


35471 


51 


1535914 


71 


34751159 


91 


529784908 


12 


272 


32 


43820 


52 


1817503 


72 


40143942 


92 


602318715 


13 


373 


33 


53963 


53 


2147434 


73 


46329631 


93 


684328892 


14 


508 


34 


66273 


54 


2533589 


74 


53419131 


94 


776998612 


15 


684 


35 


81156 


55 


2984865 


75 


61537395 


95 


881650031 


16 


915 


36 


99133 


56 


3511688 


76 


70826486 


96 


999764335 


17 


1212 


37 


120770 


57 


4125842 


77 


81446349 


97 


1132995265 


18 


1597 


38 


146785 


58 


4841062 


78 


93578513 


98 


1283193401 


19 


2087 


39 


177970 


59 


5672882 


79 


107427163 


99 


1452423276 


20 


2714 


40 


215308 


60 


6639349 


80 


123223639 


100 


1642992568 



Table 10: The number e + 1 of types of equivalence classes of FS sets of points for 
n = 1,...,100. 




sqrt(n) 



Figure 4: Plot of the data in Table |T^. 
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B Gaufi Cubature Rules for U3 



We present computed FS (usually FSMG) Gaui3 rules for Us, for m = 1, . . . , 8 (de- 
grees 3, 5, ... , 17). For m ^ 5, Mathematica (or a combination of it and some 
manual substitutions) provides analytic solutions. For higher degrees, only tran- 
scendental solutions exist. Instead, we use the matlab program cubature. m to 
approximate a solution, and we have to trust that an approximation with a small 
residual corresponds to a transcendental solution. The analytic solutions for low 
degrees provide good test data for cubature . m. All of the rules presented have been 
discovered by cubature .m; where possible, analytic solutions have been substituted. 
Beyond m = 8, no solutions at all have been found, but they should be available 
with sufficient computational effort. For most cases the rules are FSMG from the 
first structure corresponding to the first minima of the IPP. For m = 4, the struc- 
ture is from the second minima of the first structure. For m = 6 there is no FSMG 
rule, but an FSGM rule from the first minima of the second structure is obtained. 

• m = 1, degree 3. An FSMG rule is U3 : 3-1.1(1, 0, 0, 0, 0, 0)-6: 



m = 2, degree 5. An FSMG rule is C/3 : 5-1.1(1, 0, 0, 1, 0, 0)-14: 
47r 

ai = — ai = 1 

15 

Stt 1 

= To '' = W 



m = 3, degree 7. An FSMG rule is U3 : 7-1.1(1, 1, 0, 1, 0, 0)-26: 
47r 

ai = — ai = 1 

= T05 ^^ = 71 

Qtt 1 



m = 4:, degree 9. An FSMG rule is U3 : 9-1.2(1, 0, 1, 1, 0, 0)-38: 
ai = — ai = l 



97r 1 
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m = 5, degree 11. An FSMG rule is U3 : 11-1.1(1, 1,0, 1, l,0)-50: 

ai = 1 



ai = 



16n 
315 



m 



hi 


25677 

2835 ~ 


1 










277r 
~ 320 


1 








ei 


1464177 
~ 181440 ~ 


1 

\/TT 




3 

^^ = ^/n• 




6, dej 


giee 13. An FSM rule is 




13-1.1(1, l,l,l,l,0)-74: 




ai 


w 0.00644739233053 




ai 


= 1 




61 


w 0.20865289186971 




Pi 


= 1/V2 




Cl 


« 0.20762372406088 




71 


w 0.32077264898077 


Si 




« -0.37178913059595 




ei = 


= I/V3 




ei 


w 0.33396646771858 




Cl 


w 0.48038446141531 


m 


rule is not good, as di < 0. A ( 


non- 


unique) FSGM rule is U3 : 13- 


-2.1(1,0. 


ai 


w 0.05571838151106 




ai 


= 1 




Cl 


« 0.18861500631211 




7i 


« 0.33370053800545 


Ol 


ei 


w 0.12537551702973 




Cl 


w 0.70117074174860 


Vi 


62 


w 0.19567865687870 




C2 


« 0.43948383947130 


m 


7, degree 15. An FSMG rule 


is ?73 


: 15-1.1(1, 0,1,1, 2, 0)-86: 




ai 


w 0.14506632743849 




Oil 


= 1 




Cl 


w 0.14843778669299 




71 


w 0.92733065715117 


Ol 


di 


w 0.15009158815708 




ei 


= l/\/3 




ei 


w 0.13961936079093 




Cl 


« 0.36960284645415 






w 0.14924451686907 




C2 


« 0.69435400660267 




8, degree 17. An FSMG rule 


is U3 


: 17-1.1(1, 0,1,1, 3, 0)-110: 




a 


w 0.04810746585109 




a ~ 


= 1 




Cl 


w 0.12183091738552 




71 


« 0.87815891060407 


Si 


di 


w 0.12307173528176 




61 


= l/\/3 




61 


w 0.10319173408833 




Cl 


w 0.18511563534456 




62 


w 0.12058024902856 




V2 


« 0.82876998125269 


C2 


63 


w 0.12494509687253 




C3 


« 0.69042104838229 


r/3 



0.94715622136259 
0.73379938570528. 



61 « 0.94267913466612 
0.12930267526790 
0.78339511722191. 



0.37424303909034 



0.47836902881214 

0.96512403508666 

0.39568947305584 
0.21595729184587. 
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C C Code 



C.l findec.c 

/* 

Find the number of types of equivalence classes e+1 , of fully 
symmetric points in n dimensions, by exhaustively enumerating them. 
A UNIX input line "findEC a b" will find e+1 for n = a, b. 

Peter Adams and David De Wit 
August 1 1993 

*/ 

int string[100], ctr, n; 

int checkrep(ind) 

int ind ; 

{ 

int repnims [100] , rn, i, Iv; 

for (i = 0; i < 100; i++) 

repnums [i] = ; 
repnums [m = 0] = 1 ; 

Iv = 1; 

for (i = 1; i <= ind; i++) { 
if (string [i] == Iv) 

repnums [rn] ++ ; 
else { 

Iv = string [i] ; 
repnums [++rn] = 1 ; 

} 

} 

for (i = 1; i < m; i++) 

if (repnums [i] > repnums [i-1] ) 
return (0) ; 
retum(l) ; 

} 

void build (lastv, ind) 

int lastv, ind; 

{ 

int i; 

if ( ! checkrep(ind)) return; 
ctr++; 

if (ind == n) return; 

string [ind] = lastv; 
build (lastv, ind+1) ; 
if (lastv < n) { 

string [ind] = lastv + 1; 

build ( last v+1, ind+1); 

> 

string [ind] = 0; 
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main(argc, argv) 
chcir *argv[] ; 
{ 

int i; 



printfC n\t e+l\n \n" , n, ctr) ; 

for (n = atoi(argv[l] ) ; n <= atoi(argv[2] ) ; 11++) 
{ 

ctr = 0; 

for (i = 0; i < n; i++) 

string [i] = 0; 
ctr++; 

string [0] = 1; 

buildCl, 1); 

printf ("y.2d\t%6d\n", n, ctr); 

} 



C.2 ipp.c 

/* 

Exhaustively solve the IPP, using the precalculated lower bound for 
$ N $. This is a crude but fast and successful method. The data is 
output, for each m, in terms of increasing structure. Rmming the 
UNIX "sort" on the output orders it into increasing N, then structure. 



David De Wit 

August 1 — August 9 1993 

*/ 



#include <stdio.h> 



#define LIMIT 20 
#define NUMofM 21 



static int NLB [NUMofM] = 
{ 

0. 6, 14, 26, 38, 50, 74, 86, 110, 138, 162, 
190, 230, 258, 298, 342, 382, 426, 482, 526, 582 

}; 

static int c [4] [NUMofM] = 
{ 

{0, 1, 2, 3, 4, 5, 7, 8, 10, 12, 14, 16, 19, 21, 24, 27, 30, 33, 37, 40, 44}, 
{0, 0, 0, 1, 1, 2, 3, 4, 5, 7, 8, 10, 12, 14, 16, 19, 21, 24, 27, 30, 33>, 

{0, 0, 0, 0, 0, 0, 1, 1, 2, 3, 4, 5, 7, 8, 10, 12, 14, 16, 19, 21, 24}, 
{0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 2, 3, 4, 5, 7, 8, 10, 12, 14, 16} 

}; 

static double alpha [NUMofM] = 
{ 

0.00, 17.00, 7.00, 4.00, 3.00, 2.30, 1.80, 1.70, 1.45, 1.40, 1.35, 
1.25, 1.20, 1.20, 1.15, 1.13, 1.12, 1.11, 1.09, 1.08, 1.08 

}; 
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mainO 
{ 

int m, Kl, K2, K3, K4, K5, K6, N, Nvars, nstruct [NUMof M] , NUB[NUMofM] ; 

for (m = 1; m < NUMof M; m++) { 
nstruct [m] = ; 

NUB[in] = (int) (alpha [m] *NLB [m] ) ; 
for (Kl = 0; Kl <= 1; K1++) 
for (K2 = 0; K2 <= 1; K2++) 
for (K3 = 0; K3 <= LIMIT; K3++) 
for (K4 = 0; K4 <= 1; K4++) 
for (K5 = 0; K5 <= LIMIT; K5++) 
for (K6 = 0; K6 <= LIMIT; K6++) i 

N = 6*K1 + 12*K2 + 24*K3 + 8*K4 + 24*K5 + 48*K6; 
if (N > NUB[m]) 

break; 
if (NLB[m] <= N & 

Kl + K2 + 2*K3 + K4 + 2*K5 + 3*K6 >= c [0] [m] & 
K4 + 2*K5 + 3*K6 >= c [1] [m] & 
2*K3 + 3*K6 >= c [2] [m] & 

3*K6 >= c [3] [m] ) { 
Nvars = 2*K1 + 2*K2 + 3*K3 + 2*K4 + 3*K5 + 4*K6; 

printf("$ 7.2d $ & $ 7.3d $ & $ '/.Id $ & $ 7.1d $ & $ ", m, N, Kl, K2) ; 
printf ("•/.2d $ & $ 7.1d $ & $ 7.2d $ & $ ", K3, K4, K5) ; 
printf("7.2d $ & $ 7.2d $ \\\\\n" , K6, Nvars); 
nstruct [m] ++ ; 

} 

} 

printf ("\\hline\n 7.7. For m = 7.d, there axe 7.d structures in the range 7.d-7.d\n\n" , 
m, nstruct [m] , NLB [m] , NUB [m] ) ; 

> 

> 
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C.3 writestar.c 



/* 

Write down the system of moment equations as a LaTeX file. Input is 
m, and a structure K. Includes the moment data created by the 
Mathematica function USMoments.M. This code is not watertight! 

David De Wit 

August 4 — September 29 1993 

*/ 

#include <stdio.h> 

#define pi 3.141592653589793238462643383280 

main ( argc , argv) 

int argc ; 
char *argv [] ; 
{ 

int K[7], m, i, jl, j2, j3, Jl, 32, 33, eqno, pflag, nflag; 
double mom [4] [6] [11] ; 

/* Initialise some variables */ 
if (argc == 8) { 



m = 


atoi (argv [1] ) ; 








for 


(i = 1; i 


<= 6; i++) 








K[i] = atoi (argv [i+1] ) ; 






else ■ 












m = 

> 


3; K[l] = 


= 1; K[2] = 


= 1; K[3] = 0; K[4] = 


1; K[5] = 0; 


mom [O! 


[0] [0] = 


4*pi; 


mom [O! 


[0] [1] = 


4*pi/3; 


mom [0! 


[0] [2] = 


4*pi/5; 


mom [0! 


[0] [3] = 


4*pi/7; 


mom[0! 


[0] [4] = 


4*pi/9; 


mom[0! 


[0] [5] = 


4*pi/ll; 


mom [0. 


[0] [6] = 


4*pi/13; 


mom[0! 


[0] [7] = 


4*pi/15; 


mom [O! 


[0] [8] = 


4*pi/17; 


mom [O! 


[0] [9] = 


4*pi/19; 


mom [Oi 


[0][10] = 


4*pi/21; 


mom [O! 


[1] [1] = 


4*pi/15; 


mom [0! 


[1] [2] = 


4*pi/35; 


mom [0! 


[1] [3] = 


4*pi/63; 


mom[0! 


[1] [4] = 


4*pi/99; 


mom[0! 


[1] [5] = 


4*pi/143; 


mom [O! 


[1] [6] = 


4*pi/195; 


mom[0! 


[1] [7] = 


4*pi/255; 


mom [O! 


[1] [8] = 


4*pi/323; 


mom [O! 


[1] [9] = 


4*pi/399; 


mom lO'. 


[2] [2] = 


4*pi/105; 


mom [Oi 


[2] [3] = 


4*pi/231; 


mom [0! 


[2] [4] = 


4*pi/429; 


mom [0! 


[2] [5] = 


4*pi/715; 


mom [O! 


[2] [6] = 


4*pi/1105 


mom[0! 


[2] [7] = 


4*pi/1615; 


mom [O! 


[2] [8] = 


4*pi/2261 


mom [O! 


[3] [3] = 


20*pi/3003; 


mom [0! 


[3] [4] = 


4*pi/1287 


mom [0! 


[3] [5] = 


4*pi/2431; 


mom [0! 


[3] [6] = 


4*pi/4199 


mom [0! 


[3] [7] = 


4*pi/6783; 


mom [0! 


[4] [4] = 


28*pi/21879; momCO! 


[4] [5] = 


28*pi/46189; 


mom [O! 


[4] [6] = 


4*pi/12597; mom[o: 


[5] [5] = 


12*pi/46189; 


mom[l! 


[1] [1] = 


4*pi/105; 


mom[l! 


[1] [2] = 


4*pi/315; 


mom[l! 


[1] [3] = 


4*pi/693; 


mom [li 


[1] [4] = 


4*pi/1287; 


mom[l! 


[1] [5] = 


4*pi/2145 


mom[l! 


[1] [6] = 


4*pi/3315; 


mom[l! 


[1] [7] = 


4*pi/4845 


mom[l! 


[1] [8] = 


4*pi/6783; 


mom[l! 


[2] [2] = 


4*pi/1155 


mom[l. 


[2] [3] = 


4*pi/3003; 


mom[l! 


[2] [4] = 


4*pi/6435 


mom[l! 


[2] [5] = 


4*pi/12155; 


mom[li 


[2] [6] = 


4*pi/20995; mom[i; 


[2] [7] = 


4*pi/33915; 
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mom [1] [3] [3] = 

mom[l] [3] [5] = 

mom[l] [4] [4] = 

mom [2] [2] [2] = 

mom [2] [2] [4] = 

mom [2] [2] [6] = 

mom [2] [3] [4] = 

mom [2] [4] [4] = 

mom [3] [3] [4] = 



4*pi/9009; 

4*pi/46189; 

28*pi/415701; 

4*pi/5005; 

4*pi/36465; 

4*pi/146965; 

4*pi/138567; 

4*pi/415701; 

20*pi/2909907; 



mom[l] [3] [4] 
mom[l] [3] [6] 
mom[l] [4] [5] 
mom [2] [2] [3] 
mom [2] [2] [5] 
mom [2] [3] [3] 
mom [2] [3] [5] 
mom [3] [3] [3] 



= 4*pi/21879; 

= 4*pi/88179; 

= 4*pi/138567; 

= 4*pi/15015; 

= 12*pi/230945; 

= 4*pi/51051; 

= 4*pi/323323; 

= 20*pi/969969; 



/* Introduce the output */ 

printf ("Output from riuming 'writestar' : m = '/,d, ", m) ; 

printfC'K = 7.d 7.d 7.d Id 7A 7.d:\n\n", K[l] , K[2], K[3], K[4] , K[5], K[6]); 



/* Deal with Subsystem I */ 

printf ("\\be\n\7.7. Subsystem I/l:\n"); 
pflag = 0; eqno++; 
printf (" I \\[ 1 \\]\n \\eq\n "); 
if (K[l]) i 

pflag = 1; printf ("6 a_l"); 

} 

if (K[2]) { 

if (pflag) printf (" +\n "); pflag = 1; 
printf ("12 b_l") ; 

} 

if (K[3]) { 

if (pflag) printf (" +\n "); pflag = 1; 
if (K[3] > 1) 

printf("24 \\sum_-Ci=l}--[7.d} c_i", K[3]); 
else 

printf ("24 c_l"); 

} 

if (KM) { 

if (pflag) printf (" +\n "); pflag = 1; 
printf ("8 d_l"); 

} 

if (K[5]) { 

if (pflag) printf (" +\n "); pflag = 1; 
if (K[5] > 1) 

printf("24 \\sum_{i=l}-{7.d} e_i", K[5]); 
else 

printf ("24 e_l"); 

} 

if (K[6]) { 

if (pflag) printf (" +\n "); pflag = 1; 
if (K[6] > 1) 

printf ("48 \\sum_{i=l}--[7.d} f_i", K[6]); 
else 

printf ("48 f_l"); 

> 

printf ("\n \\\\\n"); 



printf ("V/.7. Subsystem 1/2 :\n"); 
for (jl = 1; jl <= m; jl++) 
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Jl = 2*jl; eqno++; 

printfC" I \\[ x-r/.d} \\]\n \\eq\n ", Jl) ; 
if (K[l]) { 

pflag = 1; printf("2 a_l Walpha.lT/.d}" , Jl) ; 

} 

if (K[2]) { 

if (pflag) printfC +\n "); pflag = 1; 
printf("8 b_l \\beta_l-{y.<i}" , Jl) ; 

y 

if (K[3]) { 

if (pflag) printfC +\ii "); pflag = 1; 
if (K[3] > 1) 

printf("8 \\sum_-[i=l>-nd} c_i \\( \\gamma_i-{7.d> + \\delta_i-{7.d} W)", K[3], Jl, Jl) ; 

else 

printfC'S c_l \\( \\ganfflia_l-{y.d> + \\delta_l~{y.d} \\) " . Jl, Jl) ; 

> 

if (KM) { 

if (pflag) printfC +\n "); pflag = 1; 
printfCS d_l \\epsilon_l-{y.d}" , Jl) ; 

} 

if (K[5]) { 

if (pflag) printfC +\n "); pflag = 1; 
if (K[5] > 1) 

printfCS \\sum_{i=l}-{y.d} e_i \\( 2 \\zeta_i-^y.d} + Weta.iT/.d} W)", K[5] , Jl, Jl) ; 
else 

printfCS e_l \\( 2 \\zeta_l-{y.d} + \\eta_l-{y.d} W)", Jl, Jl) ; 

y 

if (K[6]) { 

if (pflag) printfC +\n "); pflag = 1; 
if (K[6] > 1) { 

printfCie \\sum_-Ci=l}-{y.d} f_i ", K[6]); 

printfCWC \\theta_i-{y.d} + \\mu_i-{y.d} + Wlambda.i'ty.d} W)", Jl, Jl, Jl) ; 

} 

else 

printfCie f_l \\( \\theta_l-{y.d} + \\mu_l"{y.d} + \\lambda_l*{y.d> W) " , Jl, Jl, Jl) ; 

} 

printfC \n \\\\\n"); 

> 

/* Deal with Subsystem II */ 
if (m >= 2) { 

printf ("\y.y. Subsystem II:\n"); 

for (jl = 1; jl <= m; jl++) 

for (j2 = jl; j2 <= m-jl; j2++) { 

pflag =0; Jl = 2*jl; J2 = 2*j2; 

eqno++ ; 

printfC I \\[ x-i'm y-UAy \\]\n \\eq\n ", Jl, J2) ; 
if (K[2]) { 

pflag = 1; printf ("4 b_l \\beta_l-{y.d}" , J1+J2) ; 

y 

if (K[3]) { 

if (pflag) printfC +\n "); pflag = 1; 
if (K[3] > 1) { 

printf ("4 \\sum_{i=l}-{y.d} c_i \\( Wgamma.i'Ud} Wdelta.iT/.d} ", K[3], Jl, J2) ; 
printf C+ \\gamma_i-{y.d} Wdelta.i'-Cy.d} W)", J2, Jl) ; 
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else { 

printf("4 c_l \\( \\gaimna_l-{y.d} \\delta_l-{y.d} ", Jl, J2) ; 
printf("+ \\gainma_l-{y.d} \\delta_l-{%d} W)", J2, Jl) ; 

y 

y 

if (K[4]) { 

if (pflag) printfC "); pflag = 1; 

printfC'S d_l \\epsilon_l-{7.d}" , 31+32); 

y 

if (K[5]) { 

if (pflag) printfC +\n "); pflag = 1; 
if (K[5] > 1) { 

printfC'S \\sum_{i=l}-{y.d}- e_i \\( \\zeta_i-{y.d}- + ", K[5], 31+32); 

printf ("\\zeta_i-{y.d> \\eta_i-{y.d> + Wzeta.iT/.d} \\eta_i-r/.d> \\)", Jl, 32, 32, 31); 

y 

else { 

printfC'S e_l \\( \\zeta_l-{y.d> + ", J1+J2) ; 

printf C'\\zeta_l-{y.d} \\eta_l-{y.d} + \\zeta_l-{y.d} \\eta_l-{y.d} W)", Jl, 32, 32, 31); 

y 

y 

if (K[6]) { 

if (pflag) printfC +\n "); pflag = 1; 
if (K[6] > 1) { 

printfC' \\\\\n & & WqquadS \\sum_{i=l>-r/.d]- f_i \\( ", K[6]); 

printf ("\\theta_i-{y.d} Wnm.iT/.d} + ", Jl, 32); 

printf ("\\theta_i~r/.d> \\mu_i"{y.d> + ", 32, 31); 

printf ("\\theta_i-r/.d> \\lambda_i-{y.d> + ", Jl, J2) ; 

printf ("\\theta_i'{y.d} \\lambda_i-{y.d> +" , J2, Jl) ; 

printf ("\\iiiu_i-{y.d> \\lambda_i"{y.d} + ", Jl, J2) ; 

printf C'\\nm_i--[y.d} Wlambda.i'-Cy.d} W)", J2, Jl) ; 

y 

else { 

printfC' \\\\\n & & WqquadS f_l \\( "); 
printf ("\\theta_l*{y.d> \\mu_l'{y.d} + ", Jl, J2) ; 
printf ("\\theta_l-{y.d} \\mu_l"{y.d} + ", J2, Jl) ; 
printf C'\\theta_l-r/.d> \\lambda_l-{y.d> + ", Jl, J2) ; 
printf ("\\theta_l'{y.d} \\lambda_l-{y.d> + ", J2, Jl) ; 
printf ("\\mu_l-{y.d> \\lambda_l-{y.d> + ", Jl, J2) ; 
printf ("\\mu_l'{y.d} \\lambda_l'{y.d} W)", J2, Jl) ; 

} 

y 

printf ("\n \\\\\n"); 



/* Deal with Subsystem III ♦/ 
if (m >= 3) { 

printf ("V/.y. Subsystem III:\n"); 
for (jl = 1; jl <= m; jl++) 
for (j2 = jl; j2 <= m-jl; j2++) 
for (j3 = j2; j3 <= m-jl-j2; j3++) { 

pflag =0; Jl = 2*jl; J2 = 2*j2; J3 = 2*j3; eqno++; 

printfC I \\[ x-{y.d} y--[y.d} z'r/.d} \\]\n \\eq\n ", Jl, J2, J3) ; 

if (K[4]) { 

pflag = 1; printfC'S d_l \\epsilon_l-{y.d>" , J1+J2+J3) ; 
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y 

if (K[5]) { 

if (pflag) printfC +\n "); pflag = 1; 
if (K[5] > 1) { 

printf("8 \\sum_{i=l}-{y.d}- e_i \\( \\zeta_i-{y.d}- \\eta_i-{y.d} + ", K[5], J1+J2, J3) ; 
printf ("\\zeta_i~{y.d} \\eta_i-{y.d} + ", J1+J3, 32); 
printf ("\\zeta_i-{%d} \\eta_i-{%d} \\) " , J2+J3. Jl) ; 

} 

else { 

printf("8 e_l \\( \\zeta_l-{y.d} \\eta_l-{y.d} + ", J1+J2, J3) ; 
printf ("\\zeta_l-{y.d> \\eta_l-{y.d} + ", J1+J3, J2) ; 
printf ("\\zeta_l"{y.d} \\eta_l"{y.d} W)", J2+J3, Jl) ; 

} 

y 

if (K[6]) { 

if (pflag) printfC +\n "); pflag = 1; 
if (K[6] > 1) { 

printfC \\\\\n & & WqquadS \\sum_{i=l}-{y.d} f_i \\( ", K[6]); 
printf ("\\theta_i--Cy.d} \\mu_i~-[y.d} \\laiiibda_i-{y.d}- + ", Jl, J2, J3) ; 
printf ("\\theta_i-{y.d} \\mu_i~{y.d} \\lambda_i-{y.d} + ", Jl, J3, J2) ; 
printfC "\\theta_i-{y.d} \\mu_i--[y.d} \\lambda_i-{y.d} + ", J2, Jl, J3) ; 
printf ("\\theta_i~{y.d} \\mu_i~{y.d} \\lambda_i-{y.d} + ", J2, J3, Jl) ; 
printf ("\\theta_i-{7.d} \\mu_i-{7.d> \\lambda_i-{y.d} + ", J3, Jl, J2) ; 
printf ("\\theta_i-{7.d} \\mu_i-{7.d> \\lambda_i-^y.d} W)", J3, J2, Jl) ; 

} 

else {. 

printfC \\\\\n & & WqquadS f_l \\( "); 

printf ("\\theta_l-{y.d} \\mu_l-{y.d} \\lambda_l-{y.d} + ", Jl, J2, J3) ; 
printf ("\\theta_l-{y.d} \\mu_l-{y.d} \\lambda_l-{y.d}- + ", Jl, J3, J2) ; 
printf ("\\theta_l-{y.d> \\mu_l-{y.d} \\lambda_l-{y.d} + ", J2, Jl, J3) ; 
printf C\\theta_l-{7.d} \\mu_l-{y.d> \\lambda_l-nd} + ", J2, J3, Jl) ; 

printf ("\\theta_l--Cy.d} \\mu_l-{y.d} \\laiiibda_l-{y.d}- + ", J3, Jl, J2) ; 
printf ("\\theta_l-r/.d> \\mu_l-{y.d} \\lambda_l-{y.d> \\) " , J3, J2, Jl) ; 

} 

} 

printf C\n \\\\\n"); 
} } 

nflag = 0; 

printf ("V/.y. Also Sprach:\n"); 
if (K[l]) { 

printfC l\n \\eq\n \\alpha_l\n") ; nflag = 1; eqno++; 

} 

if (K[2]) { 

if (nflag) printfC \\\\\n"); nflag = 1; 

printfC l/\\sqrt{2}\n \\eq\n \\beta_l\n") ; eqno++; 

> 

if (K[3]) { 

if (nflag) printfC \\\\\n"); nflag = 1; 
if (K[3] > 1) 

printfC l\n \\eq\n \\gaimna_i*2 + \\delta_i'2 Wqquad i = 1, Wdots, y.d\n" , K[3]); 
else 

printfC l\n \\eq\n \\gaimna_l"2 + \\delta_l'2\n") ; 
eqno += K [3] ; 

} 
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if (KM) { 

if (nflag) printfC \\\\\ii") ; nflag = 1; 

printfC l/\\sqrt{3}\n \\eq\n \\epsilon_l\n") ; eqno++; 

} 

if (K[5]) { 

if (nflag) printfC \\\\\n"); nflag = 1; 
if (K[5] > 1) 

printfC l\n \\eq\n 2 \\zeta_i~2 + \\eta_i~2 Wqquad i = 1, Wdots, 7.d\n", K[5]); 
else 

printfC l\n \\eq\n 2 \\zeta_l~2 + \\eta_l"2\n") ; 
eqno += K [5] ; 

} 

if (K[6]) { 

if (nflag) printfC \\\\\n"); nflag = 1; 
if (K[6] > 1) { 

printfC l\n \\eq\n \\theta_i"2 + \\mu_i"2 + \\lambda_i~2 "); 
printf ("Wqquad i = 1, Wdots, "/.dXn", K[6]); 

> 

else 

printfC l\n Weq\n Wtheta_l-2 + Wmu_l-2 + Wlambda_l-2\n") ; 
eqno += K [6] ; 

} 

printf C\\ee\n\nThere are a total of 7.d equations. \n\n" , eqno); 

} 



Output from running 'writestar': m = 4, K = 10110 0: 



/[I] = 6ai + 24ci + 8di 

I[x^] = 2aial + 8ci{-f^ + Sf) + Sdiej 

I[x^] = 2aia? + 8ci(7f + ^?)+8dief 

I[x^] = 2aial + 8ci{jf + 5f) + Sdief 

I[xV] = Ac^(JfSl+JfSl)+8d^et 

I[x^y^] = Aci{jfSt+-ftSl)+8d,et 

/[x^/] = AciijfSf + jfSj) + 8diel 

I[xY] = ic,{^tSt+^iSt}+8d,el 

I[xVz^] = 8die? 

I[x^y''z^] = 8dief 

1 = ai 

1 = I'^+Sf 

1/V3 = ei 

There are a total of 14 equations. 
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D MATLAB Code 



D.l Uz Gaufi Rules 

D.l.l cubature.m 

Solve the system (*) of moment equations, for input parameters 
7. m and K, to obtain a cubature rule of degree 2m+l. 

% David De Wit 

"/, September 20 — October 27 1993 



function x = cubature(m, j, t) 
format compact; format long 



if (m == 1) 



K = 


[1 








0]; 


elseif 


(m = 




2) 




K = 


[1 





1 


0]; 


elseif 


(m = 




3) 




K = 


[1 1 





1 


0]; 


elseif 


(m = 




4) 




7. K = 


[1 





1 1 


0]; 


K = 


[1 


1 


1 


0]; 


elseif 


(m = 




5) 




K = 


[1 1 





1 1 


0]; 


elseif 


(m = 




6) 




K = 


[1 1 


1 


1 1 


0]; 


7. K = 


[1 


1 


2 


0]; 


elseif 


(m = 




7) 




K = 


[1 


1 


1 2 


0]; 


elseif 


(m = 




8) 




K = 


[1 


1 


1 3 


0]; 


elseif 


(m = 




9) 




K = 


[1 1 





1 3 


1]; 


7. K = 


[1 1 


1 


1 2 


1]; 


elseif 


(m = 




10) 




K = 


[1 1 


1 


1 3 


1]; 


7. K = 


[1 1 


2 


1 2 


1]; 



end 



M = [2 2 3 2 3 4]'; 
L(l) = 0; 
for i = 2:6 

L(i) = L(i-l) + M(i-l)*K(i-l) ; 
end 

xO = rand(K*M,l) ; 

if (K(l)), xO(L(l)+2) = 1; end 

if (K(2)), x0(L(2)+2) = l/sqrt(2); end 

if (K(4)), x0a(4)+2) = l/sqrt(3); end 

options = [1 eps eps 1] ; options (14) = 10000; 
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y. Pass the values to the solution routine. 

X = fsolveCmomenteq' , xO, options, [] , m, K, L); 



dispCm, j, t, K are:'); disp([m j t K]); 



y. Decode the resulting x into known variables: 
if (K(l)) 

a = x(l), alpha = x(2) 

end 

if (K(2)) 

b = x(L(2)+l), beta = x(L(2)+2) 
end 

if (K(3)) 

c = x(L(3)+l:L(3)+K(3)) 

gamma = x(L(3)+K(3)+l :L(3)+2*K(3) ) 

end 

if (K(4)) 

d = xa(4)+l), epsilon = xa(4)+2) 

end 

if (K(5)) 
e 

zeta 
end 

if (K(6)) 
f 

mu 
end 



x(L(5)+l:L(5)+K(5)) 
xa(5)+K(5)+l:L(5)+2*K(5)); 



x(L(6)+l:L(6)+K(6)); 



delta = x(L(3)+2*K(3)+l:L(3)+3*K(3)) 



eta = x(L(5)+2*K(5)+l:L(5)+3*K(5)) 



theta = x(L(6)+K(6)+l:L(6)+2*K(6)) 



xa(6)+2*K(6)+l:L(6)+3*K(6)); lambda = xa(6)+3*K(6)+l :L(6)+4*K(6) ) 



D.1.2 momenteq.c 
/* 

C code that compiles into a MATLAB .mex file that 'evaluates' 
the moment equations (*) . 

David De Wit and Martin Sharry and Peter Adams 
September 24 — September 29 1993 

*/ 



#include <math.h> 
#include "mex.h" 



/* Input and Output Arguments */ 



#def ine x_IN prhs [0] 

#define m_IN prhs[l] 

#define K_IN prhs [2] 

#define L_IN prhs [3] 

#def ine F_OUT plhs [0] 



static int Neq[10] = {2, 4, 7, 11, 16, 23, 31, 41, 53, 67}; 
static int M[6] = {2, 2, 3, 2, 3, 4}; 



mexFimction(nlhs, plhs, nrhs, prhs) 

int nlhs , nrhs ; 
Matrix *plhs [] , *prhs [] ; 
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{ 

double *F, *x, m, *K, *L; 
int i, mm; 

/* Check for proper nvunber of argvunents */ 
if (nrhs != 4) { 

mexErrMsgTxtC'momenteq requires four input arguments."); 
> else if (nlhs > 1) { 

mexErrMsgTxtC'momenteq requires one output argument."); 

} 

/* Assign pointers to the veirious pairameters */ 

X = mxGetPr(x_IN) ; m = mxGetScalar (m_IN) ; 
K = mxGetPr(K_IN) ; L = mxGetPr (L_IN) ; 

mm = Neq[(int) m - 1] ; 
for (i = 0; i < 6; i++) 

mm += K[i]*(M[i] - 1) ; 
F_OUT = mxCreateFulKmm, 1, REAL); 
F = mxGetPr(F_OUT) ; 

momenteq(F, x, m, K, L); 



#define pi 3.141592653589793238462643383280 
double mom [4] [6] [11] ; 

void setmomO 
i 



mom [O! 


[0] [0] = 


= 4*pi; 


mom [O! 


[0] [1] 


= 4*pi/3; 


mom [O! 


[0] [2] = 


= 4*pi/5; 


mom [O! 


[0] [3] 


= 4*pi/7; 


mom [O! 


[0] [4] = 


= 4*pi/9; 


mom [O! 


[0] [5] 


= 4*pi/ll; 


mom [Oi 


[0] [6] = 


= 4*pi/13; 


mom [Oi 


[0] [7] 


= 4*pi/15; 


mom [0! 


[0] [8] = 


= 4*pi/17; 


mom [0! 


[0] [9] 


= 4*pi/19; 


mom[0! 


[0][10] = 


= 4*pi/21; 


mom[0! 


[1] [1] 


= 4*pi/15; 


mom [O! 


[1] [2] = 


= 4*pi/35; 


mom [O! 


[1] [3] 


= 4*pi/63; 


mom [O! 


[1] [4] = 


= 4*pi/99; 


mom [O! 


[1] [5] 


= 4*pi/143; 


mom [Oi 


[1] [6] = 


= 4*pi/195; 


mom [O! 


[1] [7] 


= 4*pi/255; 


mom [0! 


[1] [8] = 


= 4*pi/323; 


mom [0! 


[1] [9] 


= 4*pi/399; 


mom[0! 


[2] [2] = 


= 4*pi/105; 


mom[0! 


[2] [3] 


= 4*pi/231; 


mom [O! 


[2] [4] = 


= 4*pi/429; 


mom [O! 


[2] [5] 


= 4*pi/715; 


mom [O! 


[2] [6] = 


= 4*pi/1105 


mom [O! 


[2] [7] 


= 4*pi/1615; 


mom [O! 


[2] [8] = 


= 4*pi/2261 


mom [O! 


[3] [3] 


= 20*pi/3003; 


mom [O! 


[3] [4] = 


= 4*pi/1287 


mom [0" 


[3] [5] 


= 4*pi/2431; 


mom[0! 


[3] [6] = 


= 4*pi/4199 


mom[0! 


[3] [7] 


= 4*pi/6783; 


mom [O! 


[4] [4] = 


= 28*pi/21879; momCO! 


[4] [5] 


= 28*pi/46189; 


mom [0! 


[4] [6] = 


= 4*pi/12597; momCO! 


[5] [5] 


= 12*pi/46189; 


mom[l! 


[1] [1] = 


= 4*pi/105; 


mom[l! 


[1] [2] 


= 4*pi/315; 


mom[l! 


[1] [3] = 


= 4*pi/693; 


mom[l! 


[1] [4] 


= 4*pi/1287; 


mom[l! 


[1] [5] = 


= 4*pi/2145 


mom[l! 


[1] [6] 


= 4*pi/3315; 


mom[l! 


[1] [7] = 


= 4*pi/4845 


mom[l! 


[1] [8] 


= 4*pi/6783; 


mom [li 


[2] [2] = 


= 4*pi/1155 


mom [li 


[2] [3] 


= 4*pi/3003; 


mom [li 


[2] [4] = 


= 4*pi/6435 


mom [li 


[2] [5] 


= 4*pi/12155; 
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mom [1] [2] [6] = 

mom[l] [3] [3] = 

mom[l] [3] [5] = 

mom[l] [4] [4] = 

mom[2][2][2] = 

mom [2] [2] [4] = 

mom [2] [2] [6] = 

mom [2] [3] [4] = 

mom[2][4][4] = 

mom[3][3][4] = 



4*pi/20995; 

4*pi/9009; 

4*pi/46189; 

28*pi/415701; 

4*pi/5005; 

4*pi/36465; 

4*pi/146965; 

4*pi/138567; 

4*pi/415701; 

20*pi/2909907; 



mom[l] [2] [7] 
mom[l] [3] [4] 
mom[l] [3] [6] 
mom[l] [4] [5] 
mom [2] [2] [3] 
mom [2] [2] [5] 
mom [2] [3] [3] 
mom [2] [3] [5] 
mom [3] [3] [3] 



= 4*pi/33915; 

= 4*pi/21879; 

= 4*pi/88179; 

= 4*pi/138567; 

= 4*pi/15015; 

= 12*pi/230945; 

= 4*pi/51051; 

= 4*pi/323323; 

= 20*pi/969969; 



momenteqCF, x, m, K, L) 
double F[], x[], m, K[], L[]; 
{ 

int iK[6], iL[6] , i, im, jl, j2, j3, p; 

double *a, *alpha, *b, *beta, *c, *gammali, *delta, *d, *epsilon, 

*e, *zeta, *eta, *f, *theta, *mu, *lambda, Jl, 32, J3; 
static int havesetmom = 0; 



if ( ! havesetmom) { 

setmomO; havesetmom = 1; 

} 

im = (int) m; 

for (i = 0; i < 6; i++) { 

iK[i] = (int) K[i] ; iL[i] = (int) L[i] ; 

} 



a = 


(double 


*) 


malice 


( 


iK [0] 


* 


sizeof (double) 


alpha 


(double 


*) 


malloc 


( 


iK[0] 


* 


sizeof (double) 


b 


(double 


*) 


malloc 


( 


iK[l] 


* 


sizeof (double) 


beta 


(double 


*) 


malloc 


( 


iK[l] 


* 


sizeof (double) 


c = 


(double 


*) 


malice 


( 


iK[2] 


* 


sizeof (double) 


gammah = 


(double 


*) 


malice 


( 


iK[2] 


* 


sizeof (double) 


delta 


(double 


*) 


malloc 


( 


iK[2] 


* 


sizeof (double) 


d 


(double 


*) 


malloc 


( 


iK[3] 


* 


sizeof (double) 


epsilon = 


(double 


*) 


malice 


( 


iK[3] 


* 


sizeof (double) 


e 


(double 


*) 


malice 


( 


iK [4] 


* 


sizeof (double) 


zeta 


(double 


*) 


malice 


( 


iK [4] 


* 


sizeof (double) 


eta 


(double 


*) 


malloc 


( 


iK[4] 


* 


sizeof (double) 


f 


(double 


*) 


malloc 


( 


iK[5] 


* 


sizeof (double) 


theta 


(double 


*) 


malice 


( 


iK [5] 


* 


sizeof (double) 


mu = 


(double 


*) 


malloc 


( 


iK[5] 


* 


sizeof (double) 


lambda = 


(double 


*) 


malloc 


( 


iK[5] 


* 


sizeof (double) 



for (i = 0; i < iK[0] ; i++) { 

a[i] =x[iL[0]+i]; alpha [i] = x[iL[0]+l*iK[0]+i] ; 

} 

for (i = 0; i < iK[l] ; i++) { 

b[i] =x[iL[l]+i]; beta[i] = x[iL[l]+l*iK[l]+i] ; 

> 

for (i = 0; i < iK[2] ; i++) { 

c[i] =x[iL[2]+i]; gammah [i] = x[iL[2]+l*iK[2]+i] ; delta [i] = x[iL[2]+2*iK[2]+i] ; 



for (i = 0; i < iK[3] ; i++) {. 
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d[i] = x[iL[3]+i]; epsilon[i] 

} 

for (i = 0; i < iK[4] ; i++) { 
e[i] = x[iL[4]+i]; zeta[i] 

> 

for (i = 0; i < iK[5] ; i++) { 
f[i] =x[iL[5]+i]; theta[i] 
lambda [i] 

} 

/* Subsystem I */ 

F[p = 0] = - mom[0] [0] [0] ; 

for (i = 0; i < iK[0]; i++) F[p] += 6*a[i] ; 

for (i = 0; i < iK[l]; F[p] += 12*b[i]; 

for (i = 0; i < iK[2] ; F[p] += 24*c[i]; 

for (i = 0; i < iK[3]; i++) F[p] += 8*d[i]; 

for (i = 0; i < iK[4]; i++) F[p] += 24*e[i]; 

for (i = 0; i < iK[5]; i++) F[p] += 48*f[i]; 

for (jl = 1; jl <= im; jl++) { 
Jl = 2.0*jl; 

F[++p] = - mom[0] [0] [jl]; 
for (i = 0; i < iK[0] ; i++) 

F[p] += 2*a[i] *pow(alpha[i] , Jl) ; 
for (i = 0; i < iK[l] ; i++) 

F[p] += 8*b[i]*pow(beta[i] , Jl) ; 
for (i = 0; i < iK[2] ; i++) 

F[p] += 8*c[i]*(pow(gainmah[i] , Jl) + pow(delta[i] , Jl) ) ; 
for (1 = 0; i < IK [3] ; i++) 

F[p] += 8*d[i]*pow(epsilon[i] ,J1) ; 
for (i = 0; i < iK[4] ; i++) 

F[p] += 8*e [i] *(2*pow(zeta[i] , Jl) + pow(eta[i] , Jl)) ; 
for (1 = 0; i < iK[5] ; i++) 

F[p] += 16*f [i]*(pow(theta[i] , Jl) + pow(mu[i] ,J1) + pow(lambda[i] , Jl) ) ; 



/* Subsystem II */ 
if (im >= 2) 

for (jl = 1; jl <= im; jl++) { 
Jl = 2.0*jl; 

for (j2 = jl; j2 <= im- j 1 ; j2++) i 
J2 = 2.0*j2; 

F[++p] = - mom[0] [jl] [j2] ; 
for (i = 0; i < iK[l] ; i++) 

F[p] += 4*b[i]*pow(beta[i] ,J1+J2) ; 
for (i = 0; i < iK[2] ; i++) 

F[p] += 4*c[i]*( 

pow(gammali[i] , Jl)*pow(delta[i] , J2) + pow(gammah[i] , J2)*pow(delta[i] , Jl)) ; 
for (i = 0; i < iK[3] ; i++) 

F[p] += 8*d[i]*pow(epsilon[i] , J1+J2) ; 
for (i = 0; i < iK[4] ; i++) 
F[p] += 8*e[i]*( 

pow(zeta[i] , J1+J2) + 
pow(zeta[i] , Jl)*pow(eta[i] , J2) + 
pow(zeta[i] , J2)*pow(eta[i] , Jl)) ; 



= x[iL[3]+l*iK[3]+i] ; 

= x[iL[4]+l*iK[4]+i] ; eta[i] = x [iL [4] +2*iK [4] +i] ; 

= x[iL[5]+l*iK[5]+i] ; mu[i] = x[iL[5]+2*iK[5]+i] ; 
= x[iL[5]+3*iK[5]+i] ; 
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for (i = 0; i < iK[5] ; 
F[p] += 8*f[i]*( 

pow(theta[i] , Jl)*pow(mu[i] , J2) + pow(theta[i] , J2)*pow(mu[i] , Jl) + 

pow(theta[i] , Jl)*pow(lambda[i] ,32) + pow(theta[i] , J2) *pow(lambda[i] ,J1) + 
pow(mu[i] , Jl)*pow(lambda[i] , J2) + pow(mu[i] ,J2)*pow (lambda [i] ,J1)) ; 

} } 

/* Subsystem III */ 

if (im >= 3) 

for (jl = 1; jl <= im; jl++) i 
Jl = 2.0*jl; 

for (j2 = jl; j2 <= im-jl; j2++) { 
J2 = 2.0*j2; 

for (j3 = j2; jS <= im-jl-j2; j3++) { 
J3 = 2.0*j3; 
P++; 

F[p] = - mom[jl] [j2] [j3] ; 
for (i = 0; i < iK[3] ; i++) 

F[p] += 8*d[i]*pow(epsilon[i] ,J1+J2+J3) ; 
for (i = 0; i < iK[4] ; i++) 
F[p] += 8*e[i]*( 

pow(zeta[i] , Jl+J2)*pow(eta[i] , J3) + 
pow(zeta[i] , Jl+J3)*pow(eta[i] , J2) + 
pow(zeta[i] , J2+J3)*pow(eta[i] , Jl)) ; 
for (i = 0; i < iK[5]; i++) 
F[p] += 8*f[i]*( 

pow(theta[i] , Jl) *pow(mu[i] , J2)*pow(lambda[i] ,J3) + 
pow(theta[i] , Jl) *pow(mu[i] , J3)*pow(lambda[i] ,J2) + 
pow(theta[i] , J2) *pow(mu[i] , Jl) *pow(lambda[i] ,J3) + 
pow(th.eta[i] , J2)*pow(mu[i] , J3)*pow(lambda[i] ,J1) + 
pow(theta[i] , J3)*pow(mu[i] ,Jl)*pow (lambda [i] ,J2) + 
pow(tlieta[i] , J3)*pow(mu[i] , J2)*pow(lambda[i] ,J1)) ; 

> } >; 

/* Also Sprach */ 

for (i = 0; i < iK[0] ; i++) 

F[++p] = alpha [i] - 1; 
for (i = 0; i < iK[l] ; i++) 

F[++p] = beta[i] - l/sqrt(2.0); 
for (i = 0; i < iK[2] ; i++) 

F[++p] = pow(gammah[i] ,2.0) + pow (delta [i] ,2.0) - 1; 
for (i = 0; i < iK[3] ; i++) 

F[++p] = epsilon[i] - l/sqrt(3.0); 
for (i = 0; i < iK[4] ; i++) 

F[++p] = 2*pow(zeta[i] ,2.0) + pow (eta [i] ,2.0) - 1; 
for (1 = 0; i < iK[5]; i++) 

F[++p] = pow(theta[i] ,2.0) + pow(mu[i] ,2.0) + pow (lambda [i] ,2.0) - 1; 

free(a); free (alpha) ; free(b); free(beta); free(c); f ree (gammah) ; free(delta); 
free(d); free(epsilon) ; free(e); free(zeta); free (eta); 
free(f); free(theta); free(mu); free(lambda) ; 



D.2 U3 Product Rules 
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D.2.1 uSprod.m 

function [y, A] = u3prod(M) 

% fimction [y. A] = u3prod(M) 
7. 

7, Product rule cubature for U_3. On input, M is the number of points in 
7. the basic rules . Gives a Gauss product rule of degree 2M - 1 on 2 M~2 
7. Cartesian points -+ y, with weights A. y is an M*2 x 3 matrix, A an 
7. M"2 column vector. To obtain a rule of degree p, a rule on 
7 (p + 1)~2 / 2 points is required. Requires gauss. m, which returns 
7 the points and weights of the $ M $-point ID Gauss-Legendre rule. 
7 See Stroud (1971), p 41. 
7 

7 David De Wit 

7 February 11 - February 12 1993 

if "existCM'), M = 5; end 

[yG, AG] = gauss(M,-l,l) ; 
j = [1:M]'; 

yCI = cos((2*j - l)*pi/(2*M)) ; ACI = ones(size(j))*pi/M; 

yl = sqrtd - yCI."2) * sqrt(l - yG.~2)'; 
y2 = yCI * sqrtd - yG.~2)'; 
y3 = ones(M,l) * yG' ; 
A = ACI * AG' ; 

y = [yl(:) y2(:) y3(:)] ; A = A(:); 

plot3(y(:,l),y(:,2).y(:.3),'+r'); 

7 Test by integrating all f mictions of x, y and z of degree < M. 
7 Works beautifully! 

7 t = M-1; m = 1; 

7 for j = : t , for k = j : t , for 1 = k : t 

7 if (j+k+1 <= t) 
7 i = [j k 1] ; 

7 tabi (m, : ) = i ; 

7 table(m) = 2* (y ( : , 1) . - (2*i(l) ) . *y ( : .2) . - (2*i (2) ) ... 

7 .*y(:,3).-(2*i(3)))'*A - ... 

7 2*prod(gamma(i+l/2))/gamma(3/2+sum(i)) ; 

7 m = m + 1; 

7 end 

7 end, end, end 
7 tabi 
7 table' 
7 m 
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D.2.2 gauss. m 

function [x, w] = gaussCn, a, b) 

% fimction [x, w] = gauss (n, a, b) 

7» Returns {x, w}, the weights and of the n-point Gauss- Legendre 
7, quadrature rule on the interval [a, b] 

Graeme Chandler 1992 

if (n==l) 

X = (a+b)/2 ; w = b-a ; return 

end 

m = l:2:2*n-l ; 

m = (l:n-l) ./ sqrt (m(l :n-l) . *m(2 :n) ) ; 
[w, x] = eig(diag(m,-l)+diag(m,l)) ; 
X = (a+b)/2 + ((b-a)/2)*diag(x); 
w = (b-a)*(w(l,:).-2)'; 
[x, m] = sort(x); w = w(m); 



y, m is off-diagonal of matrix 

y. Find spectrum of matrix 

% X are the eigenvalues 

7. w from first components 

7. ascending order 
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E Mathematica Code 

E.l IPPBasicSolution.M 

(* 

Find basic solutions to the IPP, ignoring the integer requirement, 
for m = 1, .... 20. 



David De Wit 
August 4 1993 

*) 



20 



A — 

A — 












{ 


1, 1. 


2, 


1, 


2, 


3>, 


{ 


0, 0, 


0, 


1, 


2, 


3}, 


{ 


0, 0, 


2, 


0, 


0, 


3}, 


{ 


0, 0, 


0, 


0, 


0, 


3>, 


■C- 


-1, 0, 


0, 


0, 


0, 


0}, 


-C 


0, -1, 


0, 


0, 


0, 


o>, 


{ 


0, 0, 


0, 


-1. 


0, 


0> 


} 












c = 


{6, 12, 


24, 


8, 


24, 


48> 


b = 


Transpose [ 


■c 






{ 


0, 1, 


2, 


3, 


4, 


5, 


{ 


0, 0, 


0, 


1, 


1, 


2, 


{ 


0, 0, 


0, 


0, 


0, 


0, 


{ 


0, 0, 


0, 


0, 


0, 


0, 


{- 


-1. -1. 


-1, 


-1, 


-1, 


-1, 


{- 


-1. -1. 


-1, 


-1, 


-1, 


-1, 


{- 


-1. -1. 


-1. 


-1. 


-1. 


-1, 


>] 












sol 


= Table [0, 


u. 


mmax+1}- , 



7, 8, 10, 12, 14, 16, 19, 21, 24, 27, 30, 33, 37, 40, 44}, 

3, 4, 5, 7, 8, 10, 12, 14, 16, 19, 21, 24, 27, 30, 33}, 

1, 1, 2, 3, 4, 5, 7, 8, 10, 12, 14, 16, 19, 21, 24}, 

0, 0, 0, 1, 1, 2, 3, 4, 5, 7, 8, 10, 12, 14, 16}, 

-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 

-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 

-1. -1, -1. -1, -1. -1, -1. -1, -1, -1. -1, -1. -1, -1. -1} 



nmin = Table [0, {i, mmax+1}] 
For [m = 1, m <= mmax+1, m++, 

sol[[m]] = LinearProgramming [c , A, b[[m]]] 

] 

nmin = c . Transpose [sol] 



stmp = OpenWrite["IPPBasicSolnData.tex.bak"] 
For [m = 1, m <= mmax+1, m++, 

WriteString[{"stdout", stmp}, "$ ", m-1, "$&$", nmin[[m]] 

]; 

For [n = 1, n <=6, n++, 

WriteString [{"stdout", stmp}, " $ & $ ", InputForm[sol [ [m,n] ] ] 

]; 

]; 

WriteString[ {"stdout", stmp}, " $ \\\\\n" ]; 

]; 

Close [stmp] 
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E.2 USMoments.M 



(* 

Find the moments for integrals of polynomials over U_3, and output 
into a file suitable for inclusion as a LaTeX table. The output is 
also included in the programs momenteq.c and writesteir . c . The method 
is extremely crude, but operational! 

David De Wit 

August 4 — August 12 1993 

*) 

stmp = DpenWrite ["U3MomentData.tex"] 
For [jl = 0, jl <= 10, 
For [j2 = jl, j2 <= 10-jl, j2++. 
For [j3 = j2, j3 <= 10-j2-jl, j3++, 

f[tj = (Cos[t])-(2 jl) (Sin[t])-(2 j2); 

g[t_] = (Cos[t])-(2 jl + 2 j2 + 1) (Sin[t])-(2 j3) ; 

Int = Integrate [g [t] , -ft, -Pi/2, Pi/2}] Integrate [f [t] , {t,-Pi,Pi}]; 

WriteString [ 

{"stdout", stmp}, 

"$ ", jl, " $ & $ ", j2, " $ & $ ", j3, " $ & $ ", InputFormCint] , " $ \\\\", "\n"] 

] ] ] 

Close [stmp] 
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